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ALLOCATION OF A RESOURCE 
TO ALTERNATIVE PROBABILISTIC DEMANDS: 
TRANSPORT-EQUIPMENT POOL ASSIGNMENTS* 


Rubin Saposnik and Vernon L. Smith 


Purdue University 
and 


A. R. Lindeman 


St. Louis-San Francisco Railroad 


I, INTRODUCTION 

This paper presents an analysis of a decision problem as it was found to exist on a 
large American railroad. The problem was suggested as a consequence of discussions between 
the authors and members of the management staff of the cooperating railroad concerning the 
line-haul "piggyback" transportation of truck trailers by railroad flatcar. The problem arose 
in the form of assigning flatcars and trailers to pools at the various terminals where the rail- 
road provides piggyback service, but it has application to any problem of assigning a scarce 
resource to alternative probabilistic ends. 


Il. BACKGROUND DESCRIPTION OF THE EQUIPMENT POOL PROBLEM 

Over any operating period the railroad is faced by demand requirements for outbound 
trailer loads of freight commodities to be shipped from each terminal. These requirements 
will vary, perhaps substantially, from one operating period to another, but they tend to fall into 
a general pattern, with certain terminals having relatively large outbound demands and others 
having comparatively small demand requirements. In this situation two problems arise: the 
short-run decision problem of optimally allocating the available revenue equipment to pools at 
the various terminals and the longer-run planning problem of determining the optimal stock of 
equipment to be maintained. 

In railroad parlance the expected outbound requirements at each terminal must be 
"protected" by having on hand or destined for impending arrival a "sufficient number" of 
equipment units. Hence, the concept of an equipment pool is defined as a specified stock of 
equipment units from which demand requirements are to be met. Of course, an outbound 
requirement at any terminal could be met by an inbound loaded trailer if the timing of the 
latter's arrival should happen to be such as to allow the trailer to be used to satisfy an out- 
bound order. As a matter of practice, however, this is a fortuitous event in the case of the 


*Manuscript received April 17, 1958. 
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railroad in question. Currently, the management of the railroad must "protect" 85 percent of 
its piggyback traffic by empty-trailer flatcar moves. There are two major classes of piggy- 
back service that may typically be provided by a railroad, and they have somewhat different 
effects on equipment requirements. 

Rail-billed service refers to the situation in which the railroad solicits trailer load 
lots of cargo and hires a common-carrier trucker or a rail-owned trucking subsidiary to 
perform the pickup and delivery portion of the service. In this case the railroad may simply 
rent trailers at a flat per-diem charge to provide this service. The trailers are then typically 
(though not necessarily) returned empty or, upon occasion, loaded to their point of origin. 

Truck-billed service refers to the situation in which the railroad charges a mileage fee 
for providing line-haul service for common-carrier truckers or the railroad's own trucking 
subsidiary. The major difference between the two types of service to the railroad is that under 
truck billing empty moves are revenue moves to the railroad. (As a matter of practice, rail- 
roads generally quote a somewhat lower rate to move the truckers' empty trailers.) 

In formulating an analytical statement of the pool assignment problem and in the 
empirical implementation of the model we shall work with flatcars exclusively, since applica- 
tion of the model to trailers must account for the fact that there are several different kinds of 
trailers (van, open-top, etc.) that must be pooled. Restricting the analysis to flatcars intro- 


duces considerable simplification in the application of the model without doing violence to 
principle. 


Ill. A SHORT-RUN DECISION MODEL 

Since the revenue, net of out-of-pocket expense,! which is earned in piggyback service 
tends to vary for outbound traffic from one terminal to another because of differences in the 
types of commodities moved and the average hauling distance, it would appear reasonable to 
take the maximization of net revenue as the desirable object of a flatcar assignment model.2 
For short-run applications it will be assumed that there is a fixed stock of flatcars available 
for use. The model must then assign this limited stock of equipment to pools at each terminal 
in such a way as to take into account simultaneously 
the fact that some terminals have relatively greater 
demand requirements than others and the fact that 
some terminals happen to yield higher net revenues 
per outbound trailer load. 

It will be assumed in most of the subsequent 
analysis that flatcar demand requirements at each 
terminal form a modified rectangular frequency 
distribution, as shown in Figure 1. This discrete 
density function assumes that there is an upper 
bound a; on the outbound requirements for flatcars 








at the th terminal and that the probability of having 
Figure 1 a zero flatcar requirement is p,. 





1''Out -of-pocket" or variable costs directlyassignable toa load include the pickup-and-delivery 
charges and the trailer rental. 

2In our discussions with railroad management, net-revenue criteria were used as a concrete 
means of familiarizing ourselves with the operating problem and providing management with 
some insight into our thinking processes. As indicated below, we ended with modified forms 
of the original net-revenue maximization models. 





TRANSPORT -EQUIPMENT POOL ASSIGNMENTS 


Since the minimum time required for moving a flatcar into any terminal from some 
other terminal is of the order of two days, this implies that the possible demands on the pool 
(stock) of flatcars at any terminal are compounded of the cumulative orders that may develop 
in any period of two consecutive days. Hence we let Xj denote a discrete variable representing 
the probabilistic demand for flatcars at the i th cnveidindl, where this demand is measured by 
two-day cumulations of outbound movements from the ith terminal. 

If we let Yj be the integral number of flatcars assigned to terminal i, ri be the known 
average net revenue per outbound trailer at terminal i, and S be the known fixed stock of flat- 
cars to be allocated, then the expected net profits from n terminals is given by 


n | Ji y (-p,) wi, A, (1-,) 


(1) E| n(yy--s 9) =) Zz. ities oe eee. ¢ 1. 


i=i x,=1 1 =y,+1 


The first sum in the brackets represents the expected profit from fulfilling those demands 
which do not exceed the number of flatcars assigned to the jth terminal, while the second sum 
represents the expected profit when demand exceeds the number of flatcars assigned to the jth 
terminal.3 The expression (1) is to be maximized with respect to the Yip subject to the 
constraint 


n 
(2) Z 4° 
i=1 


Equation (1) can be written in the form 
r, (1 -p;) of 


n 
E| H(yy,--+5 ¥q)| = << ». x, + yj, (a; - yj) | 5 


i= 


Since the sum of the first N integers is N(N + 1)/2, this becomes 





A, (1-p;) | 9, (9, +) 
B [risyy---1 my) |= 2 2 + y;, (a; - 


A, (1- be hale y, 1 | 
——— il, | O0<y,<a 


The relationship between expected profits at the th terminal, E (7;), and the pool assignment, 
y,; (from Eq. (1')), is shown in Figure 2. 





3The reader is referred to Chapter 9 of William Feller, An Introduction to Probability Theory 
and its Applications, or any other treatise on probability theory, for a discussion of expected 
values of random variables. 
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Figure 2 


A study of the limited data on demand requirements at the terminals of the cooperating 
railroad showed the "rectangular" distribution in Figure 1 to be a poor fit at some of the 
terminals. A better approximation at these terminals was a discrete variation on the gamma 
distribution. We define this probability (density) function over the set of nonnegative integers 
as follows: 


This last condition requires 





_(1- P)(1- e®) 
e = -B . 


e 
Hence, the probability function becomes 


Pp ifx=0 
f(x) = 


(1-P)(e?-1)e8* if x>1. 
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In most of the development of explicit solutions and in the applications to follow, we shall 
employ the modified rectangular distribution. 


IV. SOLUTION OF THE SHORT-RUN MODEL 

In maximizing (1') subject to (2) we cannot properly employ the Lagrange multiplier 
technique of the calculus, since the indivisibility of flatcars requires us to seek only integer 
values of the Vi which will yield a constrained maximum. Such integral solutions are in 
general difficult to come by, but the simple quadratic form of our objective function (1') allows 
an integral solution to be obtained quite easily for this special case. 

To demonstrate this we first substitute from (2) into (1') by eliminating y,» giving 


n-1 (1 - P;) vy, 1 
E | r(y4)--++ In| a Ne ae “2 + (ay +3) 


i=t 
% n-1 2 | 








Ss - y; | 
i n-1 
ee ew he Biles Fo 
aq | 2 a = wie 
n L \ i= | 
Necessary conditions for the existence of integers a" ssohers > eee i. , such that (3) is 


maximized, are 


oO 
E | (1°. -- Ye = er n°) | = E fi (9s. eee mn°a)| = 


(4) 
oO 
E c are % + Rs cieter Ynea)) 


BR = 3,2,.«cag R= 2: « 


That is, an integral step in any direction must not be uphill. An integer value for s is then 
determined by (2), since S must be an integer. 

If we impose the 2n - 2 conditions in (4), the following set of simultaneous linear 
inequalities is obtained: 


‘ An(L- Py) ARM -P)\ 5 r,, (1 - PL) 
b, =| —— —— Vy +-+-+ + Vy, ++--+——_—_ 9h} 


a a 
e, n 1, 





Mn (l= Py) Ag Ql ~ Py ke* & @ nye 2. 





a 
n ed 


where 
: 4, @-a,-D0G-a) 





Dy 


a 


n 
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The equality constraint (2) can, for symmetry, be written as two inequalities: 


(6) S<y,° +eeet ye ind ge” 26, 


Hence (5) and (6) together provide 2n simultaneous linear inequality restrictions on the Yi 
which in matrix notation can be stated 


gs a re) 











by; < 44, 2 eee 0 Yi s b, + ayy 
oO 
bo a Bog +++ a 0 Yo Do + ago 
(7) 
b 0 . — 
n-1 ” S ooo 8.1 a1 Yn-1 n-1 * *n-1, n-1 
s | 1 1 1 hy ~ 
where 
(1 -P,) ALCL - PD 
a = + 
an a 
and 
A, (1 - PL) 
2. 
a 
n 


This form reveals the symmetrical structure of these restrictions. 








n-1 
If we substitute S - tg = >. Yi from (2) into (5) we can obtain 
i=l 
=P) 9 Anh - Py) S- 9A) | X,(1-P,) (1 - PB) 
(8) b, < ——__—_ y, + <b + + 
% On an *y 


or 


Dy, % cue ry % (1 - PL) 
(1 _ P;,) 








@ h. a, (1 - p_) 
(s - y,°) < by, k aaa k ny, 
Ay a (1 - py) rj (1 - Py) Ae a, (1 - P,) 


(9) 
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If we sum these restrictions and substitute again from (2), then 





























n-1 a k n-l A, @ , >a) 
 - <S-y 4 Fy (8 - yy") 
kel dy (1 - B,) - Py) ~ n kel ry a (1 P;,) n 
(10) 
= -1 = 
RTM ET On (Py) 
“Lue hho) hee he OS Callie 
k=1 k k k=1 k n k 
or 
A oO A 
(11) 5 == =< Fg g«+s. 
D _ D 
where 
n-1 ),a,(S-a,-1)(1-p,) 2-1 
A = + al 
Ae @ (= 
k=1 k %n (1 - Py) k=1 
nwa J »n a, (1 P,,) 
= + 2 
a 4, a, (1- BY) 
and 
n-1 An a, (1 o P,) 
A, @ (1 - p,) 
= k on k s 
(12) w, = k=1 + : ; 
- ‘a 1 = 
e r %, (1 - P,) ; iz a, (1-p,) 
+ ee + 
— ry, (1 - p,) a iy (1 - p,) 


It follows that the integer vn must be contained in an interval of width given by Wh° 


th 


Since the terminals can be numbered in any order that we please, the n~ can be chosen so 


that a A, Ot - P,,) > ap /ry (1 - P,,) 2 0 for k = 1, 2,...,n- 1. The terminal for which the 


first flatcar yields the smallest expected profit is designated the nth, 


term in (12) will not be greater than (n - 1)/n, and the second term in (12) will not be larger 
than n - 1. It follows that 


Therefore, the first 


(13) w,<n-l. 


Hence, we obtain the important result that in general there are at most n integers ." satisfy- 
ing the condition (11). This result states in effect that in the case of multiple optima the solu- 
tions cannot involve more than n assignments to the terminal with the smallest initial expected 
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profitability. If S - A/D in (11) should happen to be an integer, then there are at most n integers 
vn satisfying (11). 


Having bounded vn in an interval containing at most n integers, we now rewrite (9) 
as follows 


Xn %_ (1 - Py) “a 
Ao (1 - P,,) 





(14) B, < y, < B, + 


oO 
. aa, G-o)G, - 4-9 
dic Hp (1 - Py) 





By = %, 


a a 


Since by construction -- > , it follows from (14) that i must be contained 
A, (1 ” P,) A, (1 - P;,) 





in an interval of width 


(15) rn %, (1 - Py) 





WwW, 4#1< 2; 


a, (1 - Py) 


Hence if the term B,. in (14) is not an integer, then there are at most two integer solutions for 


each % 


For the case n = 2, the solution above can be illustrated graphically as shown in Figure 
3. The expected profit function (1') can be mapped in the plane of yy and Yo by constructing 


Yo 
4¢ 


Iso-Expected- Profit Contours 








V4 
Figure 3 





41f B, should happen to be an integer, then there are at most three integer solutions for that 
particular y;,°. 
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iso-expected-profit contours. If we construct these contours for only integer values of yy and 
Yo» and then impose the flatcar constraint, Eq. (2), we see that a solution must lie on the 
highest iso-expected-profit contour that passes through an integral point on the constraint line, 
and there will be at most two solutions satisfying these conditions. 

The above model can be solved by an iterative method which is much superior for most 
operating applications, though it does not generate explicit integral bounds on the solution such 
as is contained in the conditions (11) to (15). The method involves the direct application of 
marginal principles by allocating successive units of the fixed stock of flatcars in such a way 
that expected profit constantly moves uphill. The first flatcar is assigned to the terminal with 
highest initial expected profit. The expected profit from adding a second car to this terminal 
is then compared with that from placing a car at each of the other terminals. The second car of 
the fixed stock is assigned to the highest of these. The expected profit from assigning the third 
car in the stock to each terminal is then compared and assigned to the highest of these. The 
process is continued until the stock is exhausted.5 Observe that this technique can be applied 
by using raw empirical distributions of demand requirements.° Applying the technique by use 
of the above modified rectangular distribution, we assign the first car to terminal k such that 
(1 - P;,) Ay 2@s- P;) div i=1, 2,...,mn. The next car is assigned to the terminal, yielding 





(1 - p;) A; ? ‘+h a. hi 
max \ i 
i 
if#k i=k 
Having assigned the first r - 1 cars in this manner to the terminals 1, 2, ...,n, Yi being the 


nonnegative number of cars assigned to the ith terminal, we assign the i~ car to terminal h, 
such that 








t 1 - p,,\ ! 1-p; : 
.-A hh  ) Ay 2 1-p; - Yj Nis Me Dy cia eg le 


V. APPLICATION OF SHORT-RUN MODEL 

As a means of testing this model, application was made to the network of eight terminal 
points between which the cooperating railroad provides piggyback service. The problem was 
to assign flatcars to pools at each of these eight terminals, assuming S = 20. 

From the data on past shipments, the distribution parameters, Pi and a; (Figure 1), 
and net revenue, Ai, were computed as shown in Table 1. 


By applying these parameter values to formulas (11) and (14) we found that optimal 
integral solutions, Ys for the pool assignments could be obtained very rapidly. These 
solution values, which happen to be unique for this case, are also arranged in Table 1. 





5ce. Richard Bellman, Dynamic Programming, Princeton University Press, Princeton, N. J., 
1957, Chapter 1, especially p. 10 and the problems on pp. 59, 108, and 139, for a discussion of 
the general problem of optimizing the sum of a set of nonlinear functions under a constraint 
such as Eq. (2). 

There may be some difficulties, however, if the marginal expected profit functions are not 
monotone decreasing. 





517742 O- 59 -2 
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TABLE 1 Concerning these optimal assignments, it 























should be observed that the average revenue at a 
Parameter terminal has a considerable influence on the size of 
Terminal the pool at that terminal. For example, terminals 3 
a, P; ri y° and 8 have the same distribution of requirements, but 
terminal 8 receives the larger pool assignment because 
1 3 | 1/2 | 117] 3 of the higher average revenue at that terminal. 
2 4) 4/5 | 78) 0 VI. A LONG-RUN PLANNING MODEL AND ITS 
3 4/ 1/2 54] 2 SOLUTION 
4 4|3/5 | 114| 3 For the purpose of determining optimal flatcar 
investment, the above model is inappropriate. The 
5 7 | 3/10) 142) 6 model to be employed for this purpose will use as its 
6 3 | 7/10} 67| 1 criterion the maximization of revenue net of variable 
7 3 | 2/5 37 | 1 costs assignable to individual loads and the investment 
cost (prorated to current account) of flatcars. 
8 4 | 1/2 | 207 | 4 If we let Y be the cost per unit time of maintain- 

















ing a flatcar (i.e., interest and repair cost) in the 
fleet, then, by use of the same notation appearing in the previous sections, expected 
long-run profit, E(7'), can be written 


2 
n d; j (1 - P,) 1 
(16) E LF (Wys-+++ Yad] = 5 a S as +3 ‘ - YY(3 


0<y, <a 


We wish to determine integral values y;° which will yield an unconstrained maximum for this 
expression. 
Proceeding as in Section IV, we must satisfy the following conditions: 


' Oo t Oo 
E fr (¥, pores Vem bec £E fr Bei iia 
(17) 
2E F" (9y°o +++ Ig # Byes Tn) 
Re FOS ey Ds 
Imposing these conditions gives 
(18) GC2£e, 2G +4; T. 


where 









~ re oa 


7a 


nm—™~™ wn ec 
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Hence, if C,. is not an integer, there is at most one integral value for ag? while if C,. is an 
integer there will be at most two integral values for » 

This explicit solution can be obtained more simply by direct application of the marginal 
conditions for profit maximization. In order for expected long-run profit from all terminals to 
be a maximum, it is necessary to choose the y; 50 that the expected profit at each terminal is 
maximized. Since the expected profit of each additional car assigned to a terminal (the mar- 
ginal profitability of that assignment) is less than that of the preceding car, the expected profit 
at any terminal will be maximized by increasing the number of cars at the terminal until the 
expected profit of the next car would be less than its marginal cost, Y. It follows that ye must 
be the largest integer satisfying 


1-p 
k Oo 
Ta, (Ok ~ %e +1), 27, 


from which the upper bound on as in (18) is obtained directly. 
It can also be shown that if demand requirements follow the modified "gamma" distribu- 
tion, the expected profit from the — flatcar at the xth terminal is 


By (y, - 1) 
(1 - P)e ki"k Ak » 


and hence x. is the largest integer satisfying 


4 = 
(1 - P,.) ok % axe 


or 


(19) y, < Z am OT 
By, Ax (1 ie P,.) 


This model was applied to the same network of eight points by TABLE 2 
use of the parameters of Table 1. Since the standard rental fee for flat- 
cars is $2.75 per day, or $5.50 per two-day period, Y = 5.5 was the value Terminal 
taken for the current account cost of flatcars. Applying the conditions 
(18) provided bounds on unique integral solutions for the %,- These 
solutions are shown in Table 2. 





tJ 
nw 
° 





Vil. THE RELEVANCE OF PROFIT-MAXIMIZING MODELS 
TO REGULATED COMMON CARRIERS 

Our co-workers at the railroad voiced considerable objection to 
the use of profit-maximization criteria for assigning railroad equipment 
to operating pools. The substance of their objection was that under ICC 
regulations 'we simply cannot pick and choose" in this way. Given two 
terminals faced with the same distribution of demand requirements, railroad operating man- 
agement does not feel free to assign more equipment to the terminal with the higher net 
revenue. 











on ouwr wn 
F&F wow sl LF LF WC WCW 
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Under ICC regulations a shipper who is faced by delays or refusals of business because 
of equipment shortages is likely to complain first to the railroad's management, but in the event 
that insufficient action is taken, he can file a legal complaint with the ICC. The outcome of 
hearings on the case may result in the railroad being ordered to purchase more flatcars. 
Analytically, such restraints would appear to modify both the short-run decision model and the 
long-run planning model. The behavior of management as it is conditioned by this regulatory 
atmosphere implies that a more reasonable criterion for short-run allocation might be the 
minimization of the expected number of flatcar shortages in the terminal system (with the hope 
of thereby minimizing the number of complaints). On the other hand, these regulations would 
seem to imply that long-run planning should be directed toward the maximization of long-run 
expected profits, subject to constraints in the form of upper probability bounds on the occur- 
rence of shortages at each terminal in the system (thereby limiting the chance of an ICC 
crackdown). 





A short-run "regulation" model will be one which minimizes E(s), the expected number 
of flatcar shortages, where 








a. 2 
= : (x, - y;) (1-p,) mn 1-p,| y; a. (a, +1) 
» “2 i i i 1\ 5 ea 
E oe = ae = : —— 
= on 2, a a a oe | 3 ( i*9/%* 3 , 
i=1 x=y,+1 i i=1 i 


subject to the constraint (2). The solution here is identical with that of Section IV, with each 
4V= | 


(21) S-—-wisy°<s-— 


where A'= A, D'=Dandw)=w, for A, =1, co? ere 


Similarly, 
(22) BL < y, < B, + W,; >. ae oe oe 


where B, = B, and W,. = Wy for 4, = 1, i= 1, 2,...,n. 


The application of this short-run regulation model to our eight-terminal system, with 
S = 20, reveals no less than four alternative optima. These alternative solutions are shown in 
Table 3. It should be noted that some of these solutions are very near to the profit maximizing 
solution of Table 1. Indeed, it would make sense for management to apply profit-maximizing 
criteria in choosing among these alternative "optima." 

This application demonstrates the importance of avoiding continuous approximations to 
problems which by nature require a discrete formulation. A continuous solution to this prob- 
lem provides a deceptively unique nonintegral solution. It is not sufficient, in principle, simply 
to round off solutions to a continuous model by some arbitrary rule which does not allow 
explicitly for the likelihood that the optimum may be a broad plateau rather than a sharp peak. 
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In view of management's objections to a straight TABLE 3 
profit-maximization model, a reasonable long-run "regu- 
lation" model might aim for the maximization of E(1') in Solution 
Eq. (16), subject, however, to the following probabilistic Terminal 
constraints onthe occurrence of shortages at any terminal.? i | @l Bel ev 
(a; te y)) (1 7 P;) 1 2 2 2 
(23) Prob {s;> o| = <M; ; 
a; 2 | 2) 8) 
3 s| Si si 3 
ie aa 4 3 2! 2! 2 
. , rr 5 5; 6] 5] 5 
where M; isthe maximum allowable probability of a shortage | 
at the th terminal. These constraints place simple lower 6 1 1 2| 1 
bounds on the pools at each terminal, that is (assuming " 2| 21 8 
M,=M,=...=M,=M), 8 3| 3| 3| 3 
( Total 
1-p,-M : 
(24) y,2 a, i ; Te Allocation | 20; 20} 20 | 20 
l- Pj (S) 























The reader may feel that imposing these probabilistic constraints does not solve the 
problem of formulating a modified profit-maximization model which allows for subjective 
managerial evaluations of equipment shortages. One might ask: 'How can you expect man- 
agement to translate these abstract subjective evaluations into probability coefficients?" 
Obviously you cannot expect this unless management knows the profit or other consequences 
of specifying different values for M. But this is precisely what the present formulation of the 
model is capable of doing! It can provide a relationship between optimum expected profits and 
the parameter M, thereby generating the necessary data to enable business policy-makers to 
appreciate the consequences of different courses of action. This parameterization of the con- 
sequences of policy alternatives enables the model to make a clearer separation between that 
portion of the planning problem which is amenable to quantitative analysis and that part which 
remains essentially an art. 

In applying this model it is necessary merely to use the lower bounds in (24) to modify 
the solution to the profit-maximization model discussed in Section VI, i.e., we satisfy the 
probabilistic constraints and then proceed to allocate cars to each terminal until the expected 
profit of the next car at each terminal would be below y. Since unconstrained long-run profit- 
maximization in our particular application gives a solution which provides a zero probability of 
shortage at all terminals, except the second, these constraints cannot modify the solution 
appreciably.8 





7See A. Charnes, W. Cooper, and G. Symonds, "Cost Horizons and Certainty Equivalents: An 
Pp y y 

Approach to Stochastic Programming of Heating Oil,"" Management Science, Volume 4, No. 3, 
April 1958, for an example of the use of probability constraints. 

By a comparison of the @; column of Table | with the values for YK in Table 2,it will be seen 
that the optimal allocation of equipment toterminals,1, 3, 4, 5, 6,7, and 8 is equal tothe upper 
bounds on the requirements at these terminals. Hence, for these terminals, the shortage 
probability is already zero under profit maximization. 
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TABLE 4 The effect of different ranges for the value of 
M on the solutions y and s°, and on optimum 
Optimal Solution, yi. expected profits E° (1'), is shown in Table 4. Itis seen 
that expected profits are not affected by specifying an 
~e 1 1 >M>0 M as low as 5 percent. Indeed, expected profit is 

20 | 20 , decreased by less than 1 percent below its maximum 
value by requiring a zero shortage probability at 
every terminal. These results are peculiar to our 
application with the modified rectangular distribution. 
The analysis, however, is easily extended to cases in 
which the modified-gamma or other discrete distribu- 
tions are used. 








Terminal 





VI. CONCLUSION 

Two basic kinds of models have been investi- 
gated, viz., those which have unconstrained profit 
maximization as their objective and those which 
involve a modification of this objective because of st 
g° 31 32 auxiliary considerations. In the context of railroad bi 
decision-making, the important auxiliary consideration F 
E° (r') 671.20 669.70 is the body of ICC regulations governing railroad com- pi 
mon carriers. The authors’ research objectives were ps 
to develop models which are, as nearly as possible, a 
analytical analogues of railroad management behavior or analytical formalizations of what ce 
management appears to be attempting to achieve. Our objectives were thus partly behavioristic 
and partly normative. The "regulation'' models of Section VII are believed to represent impor- 
tant advances toward accomplishing these objectives. 
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It has been shown that if only integral solutions to our models are considered then the 
likelihood of obtaining multiple optima becomes large. This is favorable to the possibility of 
overlap between the solutions to models with different objective criteria. Hence, management {1 
may be able, in some circumstances, to satisfy simultaneously different objective criteria which [2 
might otherwise appear to be in conflict. 

Finally, it has been shown that a profit maximization solution can be modified in the 
direction of reducing the probability of equipment shortages (providing more reliable service) 
without seriously hampering profits. This profit decline can roughly be charged as one of the 
costs to the railroad of common-carrier regulation. The model therefore provides a means of 
evaluating the costs of ICC regulation. In so far as flatcar service is concerned, our model and [4 
its application suggests that this cost is of the order of magnitude of 1 percent of net revenue. 

The context in which the general problem of allocating a scarce resource to alternative 
probabilistic demands suggested itself to the authors—the determination of optimal flatcar 
pools—is only one of many applications that might prove to be of importance. 

Other applications of such a model, in either discrete or continuous form, which might 
be made are as follows: 

1. Determination of optimal pools of any kind of equipment, e.g., boxcars, tankcars, 
trucks, taxicabs, and so forth. 


sO | 
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2. Establishment of optimal crew or driver pools at the various terminals served by an 
airline, a trucking company, a taxicab company, and so forth. These pools are often 
referred to, in transportation parlance, as "extra boards." 

3. Problems in logistics. For example, a meaningful problem faced by the Airforce is 
the construction of airplane parts packages containing an optimal number of each of 
the hundreds of airplane parts which are subject to failure. In this context the prob- 
lem is perhaps best formulated in terms of the objective of minimizing the number 
of airplane shortages due to the failure of a part not contained (or not contained in 
sufficient quantity) in the package. The constraint on choice is the maximum allow- 
able weight of the package, i.e., the sum of the weights of each of the parts times the 
number of units included in the package cannot exceed some upper bound.? 

4. Inventory problems. For example determination of optimal inventory levels in each 
of several warehouses located in the various marketing areas served by a firm. 
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AN INVENTORY POLICY FOR REPAIR PARTS” 
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1. INTRODUCTION 

This paper considers an inventory problem for a manufacturer who maintains a stock of 
repair parts of his own products. The demand for any given part is generated by failures of 
that part in the existing stock of products. Such failure distributions are known for typical 
items. In the simplest case considered here, the conditional probability of failure is inde- 
pendent of the lifetime of the part. The number of parts demanded during a time interval of 
given length is then Poisson distributed. Another simple case is that of a high probability of 
failure during the first moments of use and a constant smaller conditional probability of failure 
thereafter. In this case, the number of times a demand arises in any fixed time interval is still 
approximately Poisson, while the number of units per demand is geometrically distributed. 
Suppose that each failure leads to a demand for a repair part from the manufacturer, except in 
a constant proportion of cases (say) in which the product is scrapped. We then have a simple 
description of the demand that is exercised on the manufacturer's stock of repair parts. 

The stock of parts may be replenished by placing an order with the production depart- 
ment. Each batch of a repair part requires a certain setup cost in addition to the proportional 
cost of production. In general, processing of an order will be subject to various delays. We 
will assume that the time required for delivery is a random variable of known distribution. 
Because of this source of uncertainty, the stock department will wish to postpone specification 
of the size of the batch until the manufacturing has started, that is (roughly speaking) until the 
time of delivery. This implies that not more than one order will be outstanding at a time. 

-- Unfilled demand from customers will be backlogged at a penalty to the firm, reflecting the 
loss of good will. 

The inventory problem for repair parts as outlined is of some importance, notably in 
connection with military logistics, but it also has some methodological interest. For it falls 
into a class which is most naturally analyzed in terms of continuous time. And the study of 
inventory problems of this nature is still in its beginning. The type of stochastic inventory 
model which has been perfected in a recent monograph by Arrow, Karlin, and Scarf [1] assumes 
a division of time into discrete periods at the beginnings of which all sections are taken. "The 





*Manuscript received July 29, 1958. 
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transition to continuous time in such models will require formulation of demands as a 
continuous-time stochastic process. We believe that such a reformulation will turn out to 

be not merely more realistic but also ultimately simpler to handle analytically, but it has not 
been accomplished" [1, p. 24]. “The interested reader may refer to a recent paper on inventory 
policy for fixed delivery time by Beckmann and Muth [2] and to a monograph on queues by 
Morse [3]. 

Morse has considered the following cases on the assumption that unfilled demand is lost, 
demand is Poisson distributed, and delivery times are exponential: 

1. An order for one unit is sent whenever an item is sold [p. 139]. 

2. An order for S, the maximum stock, is sent when stock is zero [p. 146]. 

3. An order for D units is sent when the stock level falls into a reorder point s. If 

thereafter stock falls to zero, another order is placed for s such that the two out- 
standing orders equal the maximum stock § [p. 148]. 
4. Order size D is technologically fixed, maximum inventory S is an unknown multiple 
(up to 4) of this size. An order for D is sent whenever stock falls to a multiple of 
D [p. 153]. 
These policies are interesting, but not always optimal. The present paper adds an important 
case, but there remains a great number of unanalyzed alternatives. : 

In our paper on inventory policy [2], R. Muth and I have considered a simple model in 
which demand is described by a Poisson process, delivery times are constant, and order sizes 
are specified in advance. The situation differs from the present problem, because it is then 
optimal to allow several orders to be outstanding. 

To sum up our assumptions: 

1. One unit is demanded at a time. 

2. The probability distribution of demand is Poisson. 

3. The probability density of the arrival of one shipment is a function of the time 

elapsed since the placing of the order. 

4. At most, one order is outstanding. 

5. The order size is specified at the time of the delivery. 

For mathematical convenience we add the following assumptions: 

6. An optimal starting stock exists. 

7. The inventory policy is independent of the initial state of the system. 

Our first object is to derive the cost of an inventory-control system from the "inventory 
equation" as formulated for continuous time. We shall then calculate the optimal s and D for 
the case of Poisson-distributed demand, exponentially distributed lag times, and proportional 
costs of storage and depletion. 


2. THE INVENTORY EQUATION 
Notations used herein are as follows: 


x = stock 
t = time since an outstanding order was placed 
A = probability density of a demand 

u(t) = probability density of delivery 


(continuous) discount factor 
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hx for x>0O storage cost 
-gx for x< 0 shortage penalty 


k = fixed ordering cost 


S = optimal starting stock 


s = reordering point 
D=S-s 
u(x) = loss function when no order is outstanding 
v(x,t) = loss function when an order placed t units previously is outstanding 


L(s,D) = expected value of discounted cost for a starting stock s+D, when the order 
size is D and the reordering point is s 


B = @+2X 


~ 
il 


a+r +u (t) 


Suppose that stock is x and no order is outstanding. Compare the discounted expected 
cost of this state, u(x), with what will prevail after a short time interval of length At. During 
At costs f(x) - At will arise. With probability \ At, a sale will be made changing stock to 
x-1. Whether this will result in an order depends on whether u(x-1) or v(x-1,0) +k is 
smaller. Applying appropriate discount factors, we now have the identity 


-aAt -aAt 


(2.1) u(x) = f(x) At + (1-AAt)e u(x) + AAte Min [u(x - 1), v(x - 1, 0) + k]. 

Consider secondly the case where stock is x and time t has elapsed since an undelivered order 
was first placed. Then, with probability pu (t) At, the order will be delivered, changing the cost 
function to u(S); again, with probability \ At, a unit will be demanded, resulting in 


v(x-1,t+At). Thus 


aAt 


v(x, t) = f(x) At + (1-AAt - pAt)e v (x, t) 
(2.2) 
+r At set v(x-1,t+At) +pAt eat u(S) , 
In the limit, when At > 0 
(2.3) B u(x) = f(x) + >» Min [u(x - 1), v(x - 1, 0) + k], 
(2.4) Y v(x, t) eae (x, t) = f(x) +rAv(x- 1, t) + wu(S). 
ot 


Suppose that we start the system at a sufficiently high level of x and with no order 
outstanding. As x is decreased through sales, either it will never pay to order—a trivial 
case which we rule out—or a point s will be reached such that 
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Min [u(s), v(s, 0) + k] = v(s, 0) + k 
Min [u(x), v(x, 0) + k] = u(x) for all x > s .* (2. 


Equation (2.3) maynow be written more explicitly: 





Th 
(2.5) Bu(x) = f(x) + Au(x- 1), x>s+l1; 
(2.6) Bu(s +1) = f(s +1) + Av(s, 0) + Ak. anc 
Through successive substitution in (2.4), u(S) may be expressed in terms of v(s, 0) as follows: 
D=1 Sut 
1 — or. 
(2.7) u(S) =— {(S- i) a +a [k + v(s, 0)], where a = ‘ 
B é a+r 
i=o 
To calculate v(s, 0), we turn to the differential-difference equation (2.4). Multiply both sides 
by 
t t 
exp -{ (a +A +4 (ty) dty) = exp|-(a + r)t - | udt; }/=(t), say. 
0 fe) 
Suc 
Now, 
= - +(a+r( 4p) vo = -2%. re. ee 
at at at yt 
= —2 (vg). 
ot 
-@ 
(2.8) = [v (x, t) o(t)] = f(x) o(t) + u(t) o(t) u(S) + A M(t) v(x-1,t). 
Write v(x, t)¢@(t) = y(x, t) and consider this to be our new unknown. Integrating (2.8) from t We 
to © we have 
© © © 
(2.9) y(x,t) - lim y(x,t,) = f(x) @ dt, + u(S) | ugodt, +r | y(x-1,t,) dt,. 
1 1 1 1 1 
17” t t t 
For 


To estimate the limit, observe that v(x,t) cannot exceed the expected discounted cost when 
no deliveries occur. The latter is bounded by 
(2.1 





*It is a fine point of mathematics to show that this s is independent of whatever initial stocks 
x > s we started with; in other words, that there are no arbitrarily large s. For this to be 
true it is sufficient that disposal of stock is costless. For then u(x) must be a nonincreasing 
function. This implies that, from a certain stock level on, disposal actually takes place. But 
s can be no larger than this maximal stock level. Incidentally, it follows also that an optimal 
starting stock S must then exist. 
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© 
(2.10) (ts + r»te"*t " = const. + & |x]. 
a re) a 


Therefore, 


ly(x, t)| < [Cy + C4) [xl ele +t 


and 


lim y(x, ty) = 0. 
t}>0 


Substituting (2.9) with x - 1 in the last integral yields 


re) co 0 re) 
y (x, t) = £(x) | edt, + f(x-1)A J | ¢dty dt, + u(S) J uodt, 
4 t 
1 


Oo © © 0 
+ u(S) A J | uo dt, dt, +r? J | y(x - 2, ty) dts . 
t ty t ty 


Successive substitution leads to 


n : f+ 6) ie @) n : © fe 6) 
y(x,t)= | £(x-i) a! J aol odt,,,...dt, +u(s) ) ai | = dudt,...dt, 
i=o t ti ico !Clt«*t] ti 
ie8) i 8) 
ae J wets J y(x-n,t,,,) dt... --- dt,. 
n 


We now recall the formula for the iterated integral (assuming the integral to exist), 


iy Sty) dt, . 


+ 6) co 


For the last term in (2.9) we now obtain 


»2 i @) 
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By applying the estimate (2.10), it is easy to show that (2.11) is bounded by 





, 


A “a t 
(. =) nile Icy +, |x| + Cyn 


and therefore converges to zero uniformly in t(20) as n—»a. We therefore have, putting t=0, 


00 2. ti t 
y¥(x, 0) = >: f(x - i) { s" — exp| -(@ +A)t - J nat, lat 
l 0 il oO | 
i=o 
(2.12) 
cO oa) i t 
+ u(S) 2 | — +t" u(t) exp] -(a@ +A)t - { nar, |a - 
o «(il oO 
i=o 
t 
- f wet, 
Now u(t) e 7 is the probability density that delivery time equals t, which will be 


written q(t). The last term assumes, then, the form 


ioe) io 6) 


co i 
u(S) - 7. | ay qe oat q(t) dt = u(S) - J ett a(t) dt, 
izo ° : . 


after an interchange of summation and integration—which is legitimate because the sum con- 
verges uniformly in T if T denotes the upper limit of integration. 
Consider now the first integral on the right-hand side of (2.12), 


i t 
” at)’ 
——|exp -At-at - wdt, at. 
oO fe) 


il 


Now [at il ] e™= . p(i,t) is the (Poisson) probability that i units are demanded during an 
interval t. The probability that no delivery is made within a time t after placing an order is 
t 
- { u dt 
given by e ° = 1- Q(t). 
Equation (2.12) therefore takes the form 


i 8) ec 

y(s, 0) = )_f(x- i) { pay [1- ey] ent at 
: Oo 
1=0O 


(2.13) 


oO 
+ u(S) J q (t) eo ** at = v(s, 0). 
Oo 
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If the probability density for a demand depends on the time elapsed since the last demand, but 
this time is not remembered, a calculation along the same lines, but more complex in detail, 
leads to the same formula. However, p(i, t) is then no longer Poisson. 

Substituting (2.13) into (2.7) and ordering terms, we have 





D oe) (° 
1 at 
: + f(s +i) ai + {(s - i) p(i- Qe dt 
i=l 
(2.14) L(s, D) = ; sninciatiiactiigha agin caamienceaMeonabndbedae 
fe @) 
aD, | qe** at 
oO 


3. MINIMIZING COST PER TIME 
Since @ is small relative to \ and provided that q(t) is small for large t, the denomi- 
nator may be approximated as follows: 


where t + D/) is the average length of a stock cycle from §S to S. 

On the other hand, in the numerator we may suppress the discount factors, provided we 
absorb the interest cost in the charges for storage and depletion, f(x). The numerator repre- 
sents the total inventory cost per cycle. The formula, 


1 D Ca) ee) 
k +>), f(s +i) +)" (6 - i) [ pa» f- eat 
C oO 
(3.1) = ~— = L(s, D) , 


~ 00 
D 
d 





+ tq(t) dt 
Oo 


now gives the average cost of the inventory system per unit time. Up to a factor of propor- 
tionality—which is irrelevant for purposes of minimization—the formula (2.14) for the dis- 

counted cost of the system may, therefore, be approximated by the simpler expression (3.1) 
for the average cost per unit time. 
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We now specify 


x >0 


(3.2) 1 -™ 
-gx x<0O 


p(i, t) Poisson 


Writing 
oO 
t = | tq(t) dt expected delivery time 

o 

i 8) 
o* = { (t -7)* q(t) dt variance of delivery time 

) 

@ 


~ 
ul 


i J p(i, t) [1 - Q(t)] dt, 


we obtain, after some calculation, 





s-1l 
- r h h D(D + 1) >: , 
L = —|k +—sD +———- + (h + s- i) a, 
sigh r 22 ( at * 


(3.3) 
+ g— (0 + t2) - esi| . 


4. ORDER SIZE D AND REORDERING POINT s 

It remains to determine the best values for s and D. These are given by the conditions 
that the total discounted cost of the inventory system be minimal. To a first degree of approxi- 
mation the average inventory cost per unit time, L, may be minimized instead. To this end 
we set the first differences of L with respect to s and D equal to zero: 


s-l 
hD 
4, L~—-+ hee) a, - gt =0 
i=o 
s-l gt - “ 
(4.1) a, = ; 
i=o - 


To take the first difference with respect to D, write L in the form 


1 
D+a,t 


L = 





ea +hsD +2p@+| , 


A 
where k is independent of D. 
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Sin 
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Disappearance of 4, L means that 


h(s +D) (D +t) = kA +hsD +. D(D +1) 
2 


D + (At - 1/2) -) - 2srt+(rt - 1/2)? . 


Since D must be integral, terms of the order of 1/2 may be neglected. 


or 





(4.2) ae | 2KEX , (7? , 
h 
where 
‘ s-l re) 
K=k-hst=k+h )' (s-i)a, +g) (i-s)a;- hst 
1=0 1=S 
(4.3) 


8) 
r wi 
k +(e +h) DG s)ay-2X(0? 424), 
i=s 


@ 
where b (i - s) a; represents the expected size of the shortage per cycle. 
i=s 


5. POISSON DEMAND—EXPONENTIAL DELIVERY TIME 
The simplest case is, of course, that in which both \ and yw in Eq. (2.1) are independent 











s of time. This means that, during any fixed-time interval, demand is Poisson distributed and 
i- that the distribution of lag times is exponential with t = 1/u. Now 
© i i 
r ki< r 
a, = { Qty t) o* BS os : 
A O i! At+p \A + 
Aj 
= (say). 
A+pU 


Formula (4.1) for s becomes 
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(5.1) qs. Hh +4). 
g+h r 
log ~ + log ( = D) 
(5.2) gee - 
log 
At p 


To calculate D a direct approach is more convenient than the use of formula (4.2). First, 








r ft 
2 kX sp +D(D +1) 22k) w-98se4,-© . 
r »* iL 


(5.3) L= : 
ese. A+p i=o 


The difference with respect to D is proportional to 








ve : (1 - a®) 
g- 2 p29 S28 i MB eT a a 
2d 2d AN+u 2 og Uf 2 
H a 
(+n)? 


By the substitution of (5.1) and clearing terms, 


ce -— p?, (2-3) p>» g-hDt sBee, 
2X wp u? r yp? 


2Ak 1 


(D - 1/2)? = = +. 





Since D must be an integer, terms of the order of 1/2 may be neglected, 


(5.4) D= ’ 





and this is the well-known Wilson formula for optimal lot size [4, p. 33]. Notice that the order 
size turns out to be independent of the reordering point s and, in particular, of the size of the 


shortage penalty g. 


Incidentally, this formula remains valid for the more general delivery-time distribution 
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-H 
eee ee my 
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where 
La = 1 cc, > 0, 
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known sometimes as the hyper-exponential distribution. Returning to formula (4.1), we obtain 
the reordering point s in the final form 


or, explicitly, 





g+h ( J =) 
lo log {1 + a 
7 h g LL rh 


(5.5) s = - 


log (1 +4) log ( + 2) 
r r 


When lags are large relative to the spacing of demand, then yp is small compared with 
\, and we have approximately 


an oe fa OBR. 


A straightforward calculation shows that the value of the minimum i.e., the average 
inventory cost per unit time, is equal to 








L = h(s + D) +», 


or, from (5.4) and (5.6), 


alan ed er 


In the case that the failure probability of a part is high during a small initial period of 
its use and constant thereafter, we may approximate by assuming that the probability of failure 
has an exceptional value, r, in the first instant and is exponential after that. To replace one 
part will thus require n items with probability (i - r) yi The distribution of demand during 
a random interval of length t is now obtained by compounding a Poisson-distributed number of 
trials with the geometric distribution of numbers per trial. 

When both r and uw are small the solution is approximately as in the case of Poisson 
demand, with a demand density \/(1 - r). That is to say, for small probabilities of initial 
failure and long delivery times, the Poisson model applies with the appropriate demand density. 
The following case, based on data supplied by a manufacturing firm in Chicago, may serve as 
an example: 
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1 

0.01 
1.80 $ 
0.002 $ 
2.00 $ 


D Sk: = #3600 = 60 
0.002 — 


0.002 
lo 2 + log (140.01.60 
© F002 * 18 | ) 





., log 0.001 + log 1.6 
-log 1.01 





_ 7-3 +.204 _ 
-0.04 


SUMMARY 

Assume that the demand for a repair part is Poisson, that the costs of storage and 
shortage per part and of reordering are constant, that the probability distribution of delivery 
times for stock replacement is arbitrary but known, and that the order size is specified at the 
time of delivery. The inventory equation for this problem is formulated and solved; the 
optimal policy is found to be of the s,S type (Section 2). Since the loss function is closely 
proportional to the average cost per unit time (Section 3), it is found convenient to minimize 
average cost in determining order size and reordering point rather than discounted future cost 
(Section 4). The solution turns out to be simple in the case of exponential delivery time; in 
fact, the minimal order size is that given by the Wilson square-root formula (Section 5). These 
formulae are approximately valid also in the case where the initial probability of failure is 
somewhat higher than the conditional-failure probability at later times. 
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A simplified two-phase technique is presented for treating artificial 
vectors when simplex calculations are made by hand. This procedure 
allows the elimination of the extra row used to keep artificial variables 
at a zero level in Phase II, when they appear at a zero level in the basis 
at the end of Phase I. 











INTRODUCTION 

The use of artificial variables is very often necessary in the simplex method [1] in 
order to obtain easily an initial basic feasible solution. A two-phase method for solving a 
linear-programming problem when artificial vectors are needed has been presented by H. M. 
Wagner [2] and other authors. In Phase I a feasible solution is found, and in Phase II one 
searches for optimality. When some artificial vectors appear at a zero level in the basis at the 
end of Phase I (if they appear at a positive level, there is no feasible solution), it is necessary 
to annex an extra row and column to the tableau in order to maintain the artificial variables at 
a zero level in Phase II. It is the purpose of this paper to show how the addition of the extra 
equation can be avoided. 


PRELIMINARY DEFINITIONS AND OBSERVATIONS 
Assume we are trying to solve the linear-programming problem 


(1) Ad = P,,r 2 0, Py 2 0, A(mxn) , 
A = (P,,..-, P,), A= (Aq, -+-5 AQ), 
n 
(2) maximize z = >». Cy Ay 
k=1 


Suppose that artificial vectors Qi, (j=1, ..., S) associated with the artificial variables My are 


needed to obtain an mxm identity matrix for the initial basis. Denote by Api (i=1, ..., m) the 


*Manuscript received July 25, 1958. 
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variables in the current basis, and by bi the vector in the jth column of the basis matrix B. 
Any vector P, can be written as a linear combination of the basis vectors in the following way: 





m 
(3) P=) Xx %- 


Primes will be used to denote the new values of all quantities when a vector in the basis is 
changed. The rank of any matrix M will be written r(M). 
Use will be made of the following: 


1. At any stage in the simplex calculation, if Ane * 0, any vector for which Xi + 0 can 
replace b.. in the basis, even if XK < 0 and Z.- % 2 0. The new solution will be basic, since 
x, ¢ 0 and feasible since bi =p; for all i. Also z=2Z. 


2. Upon the termination of Phase I, there are three cases to be distinguished: 


a. Some artificial variables appear in the basis at a positive level. This implies 
there is no feasible solution to (1). 


b. No artificial variables appear in the basis. Then the constraints (1) are con- 
sistent and r(A) = 


c. Some artificial variables appear in the basis at a zero level. The constraints (1) 
are consistent, and there may be redundancy in the constraints (1). 


3. Suppose Phase I ends in case c. Let b. be a column of B where an artificial 
variable appears at a zero level. If some x > * 0 for one of the P,, b can be removed and 
replaced by P. while a basic feasible solution with the same z is still eetanatell (the only 
change is that Aor is now A instead of ui). If all artificial vectors can be removed in this 
way, a degenerate basic feasible solution to (1) is obtained. Therefore r(A) = m, and there is 
no redundancy. ‘ 


If the procedure described in the above paragraph leads to a condition where i. * 0 
for all P.. not in the basis and for r corresponding to the p columns of B containing artificial 
vectors, then every column of A can be expressed in terms of the m-p columns of A in the 
basis, i.e., r(A) = m-p. There is then redundancy in the constraints (1), and we know that when 
the optimality conditions are satisfied at the end of Phase II, p artificial vectors will be in the 
optimal basis. 


SIMPLIFIED TWO-PHASE TECHNIQUE 
The preceding observations allow us to simplify the two-phase method as follows: 
1. In Phase I, maximize 
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i.e., we set the prices of the legitimate variables to zero and give the artificial variables a 
price of (-1). This process can end in only one of the three ways mentioned previously.! If 
either case b or c is obtained, we go on to Phase II to determine an optimal solution. 


2. In Phase II we return to the original problem. In the last tableau of Phase I the 
zero prices assigned to legitimate vectors P. are replaced by their actual prices Ch, and the 
artificial vectors (if any remain in the basis) are assigned a price of zero. We compute the 
2, - Cy, TOW by using these new prices and use the last tableau of Phase I as the first tableau 
of Phase II. Even though some artificial variables appear at a zero level in the basis at the 
end of Phase I, no additional constraint is added to maintain the artificial variables at a zero 
level. 

Let us consider how to proceed in Phase II when some artificial variables appear in the 
basis at a zero level at the end of Phase I.* Suppose that at some iteration the artificial 
variables in the basis are at a zero level and P. is to be inserted. Take b.. to be any column 
of the basis containing an artificial vector at a zero level. 

If at least one Xk > 0, the vector to be removed by using the ordinary criterion will be 
one of the b.- In this case P, enters at a zero level, and the other artificial variables in the 
basis remain at a zero level. 

When no 0, an artificial vector will never be considered as a candidate for 
removal according to the usual criterion. If x, = 0 for r corresponding to the columns of B 
containing artificial vectors, then when some column b; of B is removed which is not an 
artificial vector, 


' Xk 
(5) br = Apr 7 > bi Xa, = App 





0, 


and the artificial variables remain at a zero level. The last remaining case is that where 
Xk < 0 when r corresponds to an artificial vector. In this situation, (5) shows that ifa 
column b; is removed which is not an artificial vector, then at the next stage the artificial 
variables can become positive. In this case we take as our rule that we remove any artificial 





lin Phasel,the same sequence of vectors will enter the basis as with Charnes' (-M) method [3] 
(except for possible differences in breaking ties). This follows because for all legitimate 
vectors 


zy, - c, = -M 2. X+, + small terms in Charnes' method; 
r 


> <2 == >= X+, + 0 in the two-phase method. 


Tr 


The r's designate the columns of B where artificial vectors occur. In both cases the vector to 
be inserted is determined by the magnitude of a Xrk- 
r 





-Theoretically, we could at the end of Phase I eliminate the artificial vectors in the way 
described previously until Xrj = 0, all j and for r corresponding to the columns of B containing 
the artificial vectors. As we shall see, these rows with Xpi = 0, all j, could then be written 
off, and we would no longer have any artificial vectors in the basis. We do not do this because 
these computations may be a waste of time in that they do not bring the present solution any 

closer to an optimal solution. 
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vector b,. with Xi, < 0 rather than b,. This gives a new basic feasible solution and keeps the 
remaining artificial variables at a zero level. 

To summarize, our algorithm avoids the addition of a supplementary constraint fo 
maintain the artificial variables at a zero level. The usual criterion is used to determine the 
vector to enter the basis. Also, the usual criterion is used to determine the vector to leave 
the basis, except in one case. If all Xk < 0 for r corresponding to artificial variables and if 
one or more x,, < 0, then one of the artificial vectors with XK < 0 is removed from the basis. 
The usual criterion would not consider the artificial vectors for removal. 

Under certain circumstances, further simplifications can be made in the computations. 
If at the end of Phase I or at any stage in Phase II, we find that XK = 0 for all k and for r 
corresponding to an artificial vector in the basis, then from (5) we see that from this point on 
the Xpk will always be zero, and the artificial vector in column r will never be removed. This 
points out that there is some redundancy in the original constraints. For this reason we can 
write off these rows as they occur and continue the iterations with a smaller-sized basis. 

Our algorithm for Phase II will lead to a different sequence of iterations than what 
would be obtained with Charnes' method when artificial vectors are in the basis at a zero level. 
With Charnes' method the artificial vectors in the basis would control the vector chosen to enter 
the basis because of the -M's. The two-phase method has the advantage in Phase TI of selecting 
the vector to enter the basis as if there were no artificial vectors present. 

Example: Let us solve the following problem: 


“Ay + Ag + 3X4 < 2 

3A, - Ag t+Ag + 2AQ < 3 

“Ay + 1.52» 7 2.529 + -SA4 = 0 
-22, + 3X» + dA. +r, =0, 
Ays Agr Agr Ag 20, 

max 3, - 2A» + Ag + Age 


The last two equations are obviously dependent. We retain both to illustrate the theory. After 
slack and artificial variables are added, the constraints become: 


~ Ay + Ag + 3AQ + Ag = 2 


Bry - Ay + Ag + 2Ag + Ag = 3 








3As a matter of fact, if one works out what the additional row does in Wagner's method, one 
finds that, when Q* should be removed (because x,” in row Q* and column P, is > 0), there 
exists at least one x. < 0 in a row corresponding to an artificial vector. The artificial 
vector can then be removed in place of Q*. Therefore, our algorithm leads to the same set of 
iterations as does the method proposed by Wagner. 
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+a, + 1.525 + 2.524 + ‘SA4 +My = 0 
- 22, + Bro + DAs +AqgtUog = oO. 


We give below the successive tableaus: 4:5 




































































c; 3 | -2 1 1 0 0 
Vectors 
fh; in Po| Pil Po | Ps | Pa | P5| Pg 
Basis 
0 Ps 2 -1 0 1 3 1 0 
= Tableau 1 
0 Ps I 3} -l 1 2 0 1 
0 Q, 0 -1 1.5 2.9 5 0 0 
0 Qo 0 -2 3 5 1 0 0 
_- - 2 -1 -1 
Zi; c; 0 3 0 0 
C; 3 | -2 1 1 0 0 
Vectors 
Cp; in Po | Py} Po | P3 | Pa | Ps} Po 
Basis 
0 Pe 3 0 3.9 8.5 3.5 0 1 
3 P, 0 1 | -1.5 -2.5|- .5 0 0 
0 Qo 0 0; 0 0 0 0 0 
z; - c; 0 0} -2.5 | -8.5]| -2.5 0 0 
cj si-=2 1 1 0 0 
Vectors 
Cp; in Po| Pi} Po | P3 | Pa | P5| Pe 
Basis 
0 P; 2.53 0/)|- .88 0 3.17) 1 | i Tableau 3 
1 P. «30 0 41 1 41) O | .12 
3 P, 88 1|- .47 0 «33; © | 2 
2; ~ c; 3 O | +1 0 +1 0 +1 






































4We do not add columns for the artificial vectors, since they are not neededafter being removed 
f from the basis. 


We are very much indebted to the reviewer for pointing out an arithmetical blunder in the 
original numerical example. He also brought to our attention Ref. [4]. 
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All artificial variables are zero initially, so that no Phase I is required. (Note that Phase I 
ends when 2" = 0; it is not necessary that the optimality conditions for Phase I be satisfied.) 
At the first iteration our special rule must be used, because the usual rule would remove Pe 
and the artificial vectors would become positive. It was arbitrarily decided to remove Q;- At 
the second stage we can cross out the Qo row, since all Xe = 0. 
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ON QUADRATIC PROGRAMMING* 


E. M. L. Beale 
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The algorithm for quadratic programming given by Beale [1] is 
described indetail. Some extensions of the algorithm are indicated, and 
the relationship between this algorithm and that given by Wolfe [5] is 
discussed briefly. 











1. INTRODUCTION 

This paper is concerned with an algorithm given by Beale [1] for minimizing a convex! 
quadratic function of variables subject to linear inequalities. This algorithm does not seem to 
have attracted too much attention, perhaps because it was written up over-concisely. A more 
detailed exposition of the algorithm is given in Section 2 below. It is hoped that this will be of 
interest, both as one efficient method, though not the only one, of solving a problem of some 
practical importance, and also as an illustration of the power of Dantzig's simplex method. 

Some extensions of the algorithm are considered in Section 3. In Section 4 some 
remarks are made about the relationship between this algorithm and what seems to be its 
chief competitor for this problem, an algorithm given by Wolfe [5], which also uses the simplex 
method. 


2. THE ALGORITHM 

To describe my algorithm, I will start by outlining the simplex method in terms suf- 
ficiently general to cover my application to quadratic programming. 

Suppose that we wish to minimize some convex function C of m variables x, that must 
be nonnegative and satisfy the m (<n) linearly independent equations, or constraints, 


2D, Ays = By: OS @ bh seve Qs 


n 
j=1 





*Manuscript received November 28, 1958. 

Research sponsored by the Office of Ordnance Research, U. S. Army, Statistical Techniques 
Research Group, Princeton University, where this paper was prepared. Dept. of Army Proj. 
5B 99-01-004; Ord. R&D Proj. PB2-0001; OOR Contract 1715, Army Contract DA 36-034- 
ORD- 2297. 


1A convex function is one that is never underestimated by linear interpolation,i.e., is such that 
f£Ax + (1 -A)y) SZ AL(x) + (1 -A) fly) if OLA SLI. 
Actually, the assumption of convexity is used only to ensure that a stationary value of C is 


necessarily an absolute minimum. This point is discussed at the beginning of Section 3. 
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If these constraints are not inconsistent, we can find a basic feasible solution in which 
m of the variables take positive values¢ and the others are zero. We can then use the con- 
straints to express the basic variables, i.e., the ones that are positive, in terms of the non- 
basic variables. We will then have 


n-m 
m= Mag + 2, My Be (h = jys-e+r i) 
{21 
where z, denotes x. ‘ 
4 Jm+2 
It is customary to write these equations in tableau form, with the z» on the left-hand 
side of the equations, so that the ays are recorded with opposite signs; but I find it easier to N 


think of the problem in this form. I call the set of variables on the right-hand sides of these 
equations the nonbasic set. 
In the present trial solution of the problem, the basic variables Xp, equal ano’ and the 
nonbasic variables equal zero. 
We can now use the above equations to express the objective function C in terms of the v 
nonbasic variables and can consider the partial derivative of C with respect to any one of them, 
say Z,, assuming that all the other nonbasic variables remain fixed and equal to zero., 
If now 3C/d 2 > 0, then a small increase in Zi» with the other nonbasic variables held 
equal to zero, will not reduce C. But if 3C/a Z, < 0, then a small increase in 2, will reduce 
C. If C is a linear function, then 9 C/2 2, is constant, and it will therefore be profitable to go 
on increasing 2, until one has to stop to avoid making one of the basic variables negative. a 
If C is a nonlinear function with continuous first derivatives, then it will be profitable iz 
to go on increasing 2, until either 


(a) one has to stop to avoid making some basic variable, say X negative, or si 


(b) 3C/a 2, vanishes and is about to become positive. 


Vi 
In Case (a), which is the only possible case when C is linear, one changes the basis by tk 
making x, a nonbasic variable in place of Ze» and uses the equation tl 
r 
n-m 

i = We * Fa ye Fe r 
Lal ke 
ky 
to substitute for z, in terms of xX, and the other nonbasic variables throughout the constraints re 
and also the expression for C. al 

In Case (b), we are in trouble if C is an arbitrary function. But if C is a quadratic 
function, then 2 C/a Zy is a linear function of the nonbasic variables. We can then represent C h: 
by a symmetric matrix (cy») for k, 4 = 0,1, ...,n-m, such that p) 








2We may have to use €-perturbations to avoid degeneracy in the usual way (see Charnes [3]). 
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where Z, = 1 and z,,...,2Z,_,, denote the nonbasic variables. 


This means that 


n-m n-m n-m 
Ce. +87.  % + 7. 7. Cy p By Zp - 
k=1 k=1 f=1 


Now 
n-m 
. 2. c + C) pZ 
9 = “ko } } kt “2° 
2 dz, | 


We can therefore introduce the new nonbasic variable 


n-m 


and use this equation to substitute for 2 throughout the constraints and the objective function 
in terms of u, and the other nonbasic variables. 

Note that if 2, is an x-variable we will have an extra basic x-variable after this 
substitution. 

The only slightly peculiar thing about this is that Uy, is not restricted to nonnegative 
values. I therefore call it a free variable, calling the x-variables restricted variables. But 
there is no objection to having free variables in the nonbasic set. One simply has to remember 
that if 2 C/ dup. > 0 then C can be reduced by making u, negative—or, alternatively, by 
replacing U, by the variable Vv, = -t, and increasing Yi in the usual way. 

There is one other point about these free variables. Once a free variable has been 
removed from the nonbasic set one can forget about it. There are only two reasons for 
keeping track of the expressions for the basic variables in terms of the others. One is to 
know their values in the final solution when it is obtained. The other, and more fundamental, 
reason is to prevent them from surrepticiously becoming negative. Neither reason applies to 
any basic free variable in this problem. 

There remains the problem of motivating this particular choice of free variable. We 
have seen that it is convenient to have a set of nonbasic variables that are all zero at the 
present trial solution. But any expression of the form 


3 Og Ss 


c a. a. a 
ko kk “k tek 


would satisfy this condition, and it might seem much simpler to put all a, = 0. 
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I think the real answer is that the transformation actually made would be the natural 
one if one were minimizing a convex second-degree polynomial with no inequality constraints. 


The point is that when we have substituted Uy, for 2 inthe expression for C we willhave aterm 
ont ue plus a convex second-degree polynomial in the remaining n - m - 1 variables. The 
variable uy, only appears in this one term. So we have effectively reduced the dimension of 
the problem by 1, subject to the condition that uy, = 0. 

In the absence of inequality constraints, every change of basis would introduce a new 
nonbasic free variable, and at the end of n - m steps we would have C expressed in the form 

n-m 
Coo + >: Chi ue - The minimum value of C is then Coo? taken when all the Uy, vanish. 

k=1 

When there are inequality constraints we cannot guarantee that the algorithm will 
terminate in n - m steps, though I will show later (following [1]) that it must terminate in a 
finite number of steps. But first I will illustrate the algorithm with a small-scale numerical 
example. This example has been constructed to display all the features of the algorithm, 
though this involves giving the misleading impression that the algorithm is normally rather 
slow, even on such small problems. 

Consider the problem of minimizing 


C = 6 - 6x, + 2x? - 2x, xo + 2xz, 
subject to the constraints 
x, 20, X92 0, X, + Xq < 2; 
i.e., KX, 20, X,20, x,20, xX, + Xo + XQ = 2. 


Let us assume that we start the algorithm from the trial solution Xy = Xq = 0. But 
first let us consider the geometrical interpretation of the progress of the algorithm (Figure 1). 

Starting from the point a, we first increase X1> keeping Xo = 0, until C starts to 
increase again, i.e., until we reach a horizontal tangent to some curve C = constant. We are 
then at the point b, and introduce the new nonbasic variable Uy, defined to vanish at all points 
where a curve of constant C has a horizontal tangent, i.e., all points where C cannot be 
reduced, if Xo is kept constant. 

We now find it profitable to increase Xo, keeping ay constant. In the absence of the 
constraint Xg 20, we could have increased Xo until we reached the unconditional minimum of 
C (at X, = 2, Xo = 1), and this would be the end. "Unfortunately" we are stopped at the point c. 
This is not the final solution, because we are no longer interested in points where curves of 
constant C have horizontal tangents, since the boundary we are now on slopes down at an angle 
of 45°. We therefore move from c to d, replacing the nonbasic free variable uy by the non- 
basic free variable Uo; defined to vanish at all points where C cannot be reduced, keeping Xq 
constant. And this is the final solution. 

Let us now follow through the algorithm numerically. 

We wish to make Xy and Xo the nonbasic variables, so we write our constraint as 


Xg = 2 - Xy - XQ. 
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Figure 1 








We must express C in terms of the nonbasic variables. This is already done. The 
matrix (Cc) ») is 


6 -3 O 
4 2 «1 
> = 2; 


but I prefer to write this in longhand as follows: 
C=(6- 3X, ) 
+ (-3 + 2x, - Xo) xy 


+( - X, + 2X5) Xo. 


1 
Now the partial derivative of 2° with respect to Xy is -3. It is therefore profitable to 


increase X14) and this remains so until 
-3 + 2x; - X» = 0, 


i.e., until Xy 


372. It remains feasible to increase x, until 
ye = 3 - X - XQ 20, 


i.e., until xX, = 
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We therefore continue until X, = 3/ 2, and introduce the new nonbasic free variable 
uy = -3 + 2x, - Xp 


in place of Xy- 


We must now substitute for Xy in terms of Uy and the other nonbasic variables. It is 
convenient to do this in two stages. We first substitute for the Xy in the constraints and within 
the rows of the expression for C; and then substitute for the Xj multiplying one of the rows in 
the expression for C. We therefore have 


fo 


» 2s i 

x, = -*em * o%2: 
—— 

Xa = _* u- Xo » 
" 3 

c+ (dad 
+ ( Uy ) x, 


po loo 
1 


oY + 5%) * 
“ty Tel de 


" 
- 
nw | 
1 
bo loo 
Pad 
ee 


“_ 3x.) 
+ 9 +72 2° 


It is a straightforward matter to verify that the final expression for C will be symmetrical if 
the original one was. 

We are now at the point b in Figure 1, but before continuing I want to emphasize how 
similar the change of basis operation is to the normal change of basis operation in linear 
programming by the simplex method. We choose the variable to be removed from the nonbasic 
set by finding one associated with a negative coefficient in the first row (or column) of the 
expression for C. The choice of variable to replace it in the nonbasic set is exactly the same 
as in linear programming, except that we have another possible candidate—the negative of the 
row of C, multiplying the variable being removed. 

Finally, the pivotal operation starts out just the same. We operate on the constraints 
and then on the rows in the expression for C. Only at the end is there a difference, when we 
have to continue by operating on the columns of the first transform of C in the same way. 













i.e., until Xo = 
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Returning to the example, we find it profitable to increase Xo until 


1. Unfortunately, this is not feasible because Xq < 0 if Xo 


fore make Xq nonbasic in place of Xo, writing 


~ 


, 2 
: $4 
- S444, 
; 3 
a i 
= - +s% 
1 
+ ( -"t 
~ (-1 ~ -" 
«fias 
$+ 8s? 
+ ao 
+ (2, 3, 
°° 2 


Xo ; 


Xa 


w |p 
* 
wo 
Sa 


win 
Mad 
wo 
neal 
Mal 
io) 


We are now at the point c, and must continue by decreasing uy until either 


Us 
or Xo 
or xy 


The first of these events is 


have 





1 
3 


«2 
3 


Us = 0, so we make Us nonbasic in place of uy: We then 


2 
+a * X3 = 0, 
1 
ant, te X, = 0, 
1 
+s *3a°* 


1 3. We there- 
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“en i 
4% <2 + oe *S%: 
—— i 
“2 9° 9027 9%3” 
3. ne 1 
x, = .*s. @ 


C= (= + Say += 
2 2 2 


Ei + 
9 2 


| 
fal 
w 
sagen 


+ 
, 
bp le 
+ 
nw le 


iT 
== Fo 
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3 
2 Us Us 


(: 
+ = + 
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It is now not profitable to change either Uy OF Xa, and since C is convex and dif- 


ferentiable, no change in Us and Xg together can be profitable, and we have the optimal 
solution, given by 


np le 
Pad 
wo 
i 
Pad 
wo 


= * S, =” 
2 

Note that we needed the expression for Uy to substitute for Uy in the previous tableau; but if 
we had to make another change of basis we would ignore the expression for uy in all further 
work, 

Let me now prove that the algorithm must terminate in a finite number of steps. 
Actually, I can do this only if I make the rules of procedure a little more specific. I insist 
that if there are any nonbasic free variables u for which 0 C/du # 0, then I must remove a 
free variable from the nonbasic set. 

I also need a further definition. I say that the objective function C is in standard form 
if the linear term contains no free variable. When C is in standard form, its value in the 
present trial solution, Coo is a stationary value of C, subject to the restriction that all the 
present nonbasic restricted variables take the value zero (without any restrictions on the signs 
of the basic variables). So there can be only one such value for any set of nonbasic restricted 
variables. Now, if we make the usual €-perturbations to avoid degeneracy (applied to the 
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original constraints and not to the objective function, so that the free variables are unaffected), 
we know that C decreases at every step. So it can never return to a standard form with the 
same set of nonbasic restricted variables, even with a different set of nonbasic free variables. 
There is only a finite number of possible sets of nonbasic restricted variables, so the iteration 
must terminate if it always reaches a standard form in a finite number of steps whenever it is 
not already in standard form. I will now prove this. 

Our rules of procedure ensure that when C is not in standard form a free variable will 
be removed from the nonbasic set. So s, the number of nonbasic free variables, cannot 
increase. Further, if the new nonbasic variable is free, it is easy to show that the off-diagonal 
elements in the new expression for C in the row and column associated with the new nonbasic 
variable must vanish. It follows that C does not then contain a linear term in this variable, 
and furthermore C can never contain a linear term in this variable unless some other 
restricted variable becomes nonbasic, thereby decreasing s. 

Therefore, if C is not in standard form, and s = Sy» Say, then s cannot increase, and 
must decrease after at most So steps, unless C meanwhile achieves standard form. Since C 
is always in standard form when s = 0, the required result follows. 


3. EXTENSIONS 

It is of interest to ask what the above algorithm can do if C is quadratic but is not 
known to be convex. It is clear that it may produce a local minimum that is not a global 
minimum, and in fact it could conceivably produce a stationary value of C that is not even 
a local minimum. Two further rules improve the algorithm in this respect. One is that a 
free variable should be removed from the nonbasic set if C contains any off-diagonal terms 
in this variable, even if the linear term happens to be zero. The other is that a free variable 
should be removed from the nonbasic set if C contains no off-diagonal term in this variable, 
but the diagonal term is nonpositive—in this case the variable may be either increased or 
decreased until some basic restricted variable becomes equal to zero. With these additional 
rules, the algorithm must produce a genuine local minimum unless the linear term in some 
restricted variable in the final expression for C happens to vanish. In this case the final 
solution must be a virtual local minimum, by which I mean a point that could be made into a 
local minimum by an arbitrarily small change in the coefficients of the objective function. 

Recently Dr. Wolfe suggested to me that practical quadratic programming problems 
often involve many variables with only linear terms and a much smaller number of variables 
with quadratic terms in the objective function. This means that the rank r of the matrix of 
the purely quadratic part of the objective function might be much less than its maximum pos- 
sible value, n - m. Now this rank will be unaffected by any change of basis. And this proves 
that (if the two extra rules described above are incorporated in the algorithm) one can never 
have more than T nonbasic free variables at any one time. For before one could introduce an 
(r + 1)8t such variable, C would have to be in standard form with r nonbasic free variables. 
There would therefore be r nonzero diagonal elements in the quadratic part of C with no 
nonzero off-diagonal elements in the same row or column. This implies that all the remaining 
elements of C referring to quadratic terms must vanish if the rank is to be not more than r. 
And this in turn implies that the next new nonbasic variable (if any) must be a restricted 
variable that was previously basic. 

Consequently the total number of basic restricted variables cannot exceed m+r at any 
one time, and hence the number of rows in any tableau cannot exceed m+r, plus 1 fora 
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possible outgoing nonbasic free variable, plus n - m + 1 to record the elements of C, making 
n+r +2 in all—as opposed to the 2n - m + 2 that may be required in general. 

Indeed, if r is much less than n - m, it may be profitable to express C, not bya 
square matrix of order (n - m +1) but by a linear term 2X5) plus r squares of linear 
functions Ay ees A. . Hence 





and we record do Ay a r.. as linear functions of the nonbasic variables. If C is not convex, 
then to avoid complex numbers we may have to write 


This would reduce the maximum number of rows in any tableau to m + 2r + 2. It would 
also simplify the pivotal operations, since the change of variables in the expressions for 
Ao» > pro r , would follow precisely the same pattern as in the expressions for the basic 
variables. The one feature of this formulation that is unfavorable if r is at all large is the 
computation of the expression for =" C/ OZ, . This then involves multiplying each i by the 


coefficient of Zs, in this expression, summing over i=1,..., r, and adding in the coefficient 
of 2 in ro to the constant term. 


The next topic is an application of quadratic programming rather than an extension, but 
since it involves some special features I will illustrate it briefly. 

Much attention has recently been paid to problems of finding a minimum, or at least a 
local minimum, of an arbitrary objective function of several variables, with or without inequal- 
ity constraints on the possible values of the independent variables. The most popular approach 
seems to be to fit a linear or quadratic approximation to the objective function in the neighbor- 
hood of the present trial solution and to use this to extrapolate for a better trial solution. 
Another approach is to consider some cuboid (i.e., hyper-rectangle) in the space of the 
independent variables that is believed to contain the required minimum of the objective function. 
A quadratic approximation to the objective function can be fitted to all previously computed 
values of the objective function at points within the cuboid; and the objective function can then 
be evaluated at the point where the quadratic approximation is minimized within the cuboid. 
This evaluation will then indicate whether the cuboid should be shifted, or contracted, or both 
before the process is repeated. Such an approach can be expected to reduce the number of 
function evaluations required to give an adequate estimate of the minimum, at the expense of 
considerable extra computation between evaluations. It may therefore be appropriate if each 
evaluation is a costly or time-consuming operation. The details of such a scheme are con- 
sidered by Beale [2]. Here I am only concerned with the problem of minimizing the quadratic 
(approximate) objective function over a cuboid. The special features are the usual bounded 
variables technique and the fact that it is convenient to take our first trial solution at some 
point inside the cuboid, since the objective function will normally be computed in terms of 
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So our original nonbasic variables will be free variables. 
A simple example is illustrated below. 


Minimize 
C = 10 - 4x) ~ 2X5 + 4x? - 2x, x 
subject to 
- 25 x,5 + 1,-1< x9<+ 
Writing 
xy = xX, + 2, Xo =X) + 1, 
we have 
0 xy = 2 + X 
0 Xo = 1 + Xp 
Cc = ( 10 ad 2x1 = Xo) 
+ (-2 + 4X, - Xp») Xy 
+ (-1 - Xy + Xp) Xo 
Hence, if 
= * - 2 e 4x, - Xp, 





we have 





coordinates with their origin at the previous minimum computed value of the objective function. 


oe, 
z. 
$3, 
<3, 
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We now want to increase Xo, and the first obstacle occurs when Xo = 1, since Xo 
















reaches its maximum permissible value. We therefore write 


x5 = 2 - (1 + Xo) ’ 
i.e., 
x = 1 ¥ x, ’ 
2 2 hai 
- sol 
- 11 1 1 + 
O<x, = — 4+ —-U, - —X <3, for 
1 4 4 1 2 
O< Xo = 2 - Xo <2, 
C = (= mm = x) | 
2 2 ori 
gr 
- ju 
i ( 4 1 1 If t 
we 
3 + wh 
+ [-= ~ =i 
( 4 2) ail Th 
te1 
ba: 
“ (2 Fa +x) 
4 
re: 
1 pre 
- ( Pe | Ju, de: 
is 
int 
«2 + Sah. 
4 4 


This is the final solution, C = 27/4, where xy = 11/4, Xo = 2. We must now refer back 
to the original problem to find x, and x» in terms of xX} and Xo. In fact, x, = xy - 2= 3/4, 
Xo = Xo » 1 = a 


I now consider what to do if, after solving some quadratic programming problem, one 
wished to modify the problem by adding a further linear inequality constraint on the variables. 
This is a very simple matter when the objective function is linear, and is in effect a basic 
operation in the dual simplex method for solving linear problems. In quadratic programming 
this operation seems to be much more difficult, but a general method of performing it will be 
illustrated by means of a detailed commentary on the process of adding the further constraint 
Xx £1 after reaching the final solution to our first example. 

Our first task is to use the expressions for the basic variables to express our new 
constraint in terms of the nonbasic variables. We therefore have 
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X4 = 1 - x,, t 
— ae 
2 22 923° 


If the constant term in this expression had been nonnegative, the new constraint would 
have been satisfied at our present trial solution, which would therefore have been the required 
solution. If this is not the case, we add a parameter t to the right-hand side of the expression 
for the new constraint: 


Our present trial solution is optimal and feasible for t equal to the negative of the 
original constant term in the expression, i.e., t = 1/2, and we wish to solve the problem for 
gradually diminishing values of t until we reach t = 0. 

For t just less than 1/2, X4 < 0, and to remedy this we make X4 a nonbasic variable. 
If the expression for the new nonbasic variable contains any nonzero term in any free variable, 
we make such a variable basic. If there are several such free variables it does not matter 
which is chosen, since they must all be removed from the nonbasic set before we can proceed. 
These additional free variables will all be replaced by other free variables, since the linear 
term will be arbitrarily small for this value of t, so there is no danger of making another 
basic restricted variable negative at this stage. 

Special mention should be made of the degenerate situation where more than one basic 
restricted variable vanishes as it is about to go negative for the same value of t. A safe 
procedure is to use the e-perturbations to break the tie and return C to standard form after 
dealing with one of these variables before considering the next. But the number of iterations 
is likely to be reduced if as many as possible of these variables are made nonbasic without 
intermediate substitutions of new nonbasic free variables for old. 

In our example we have 
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6t 


(- 1 3t 


We now have C in standard form, so we can consider the effect of reducing t below 
its present value of (1/2) -¢. We must stop if the value of either some basic restricted 
variable or some partial derivative of C drops to zero. Here we are held up by the fact that 
(1/2) dC/ AX, = 0 when t = 1/3. We therefore put 


and we have 
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We can now reduce t to zero without any basic restricted variable or any partial 
derivative of C becoming negative, so we have the final solution 


c=3 
2 


3 = = 
when x, = 1, Xo = de ley 


1 
2 

It is clear that the same procedure will solve a problem considered by Wolfe [5], where 
a linear function multiplied by a parameter t is added to the objective function, and one 
requires the optimal solutions for every t. These can be traced out by gradually increasing 
(or decreasing) t and changing the set of nonbasic variables whenever some basic restricted 
variable or some partial derivative of C becomes negative. The only unsatisfactory feature 
is that when some new nonbasic restricted variable is introduced, we have to reduce C to 
standard form before we can proceed—and this in general involves replacing all the present 
nonbasic free variables. 

Note that the parameter t enters linearly into the constant terms in the constraints and 
in the linear terms in the expression for C. It also enters quadratically into Coo? but is absent 
from all other terms. 

In principle, the same approach can be used to solve a general parametric quadratic 
programming problem, i.e., with quadratic factors in C multiplied by t. The coefficients are 
then, in general, rational functions of t, which are unpleasant, though not impossible, to handle. 


4, WOLFE'S METHOD 

The algorithm described in Section 2 is very similar in effect to that presented by 
Wolfe [5]. Wolfe's method is based on the Kuhn-Tucker conditions for a stationary value of a 
function subject to inequality restrictions, developed in [4]. A heuristic derivation of the "short 
form" of Wolfe's method can be given as follows. 

By using matrix notation, let the constraints be 


Ax=b, 
and let the objective function be 


2px+x Cx. 


Then one's first thought might be that a necessary condition for a minimum is that the 
partial derivative of the objective function with respect to each x-variable must vanish, i.e., 


p+Cx=0. 


But one's next thought should be that this is wrong, because of the constraints. So one 
introduces Lagrange multipliers, and has the condition 


p+Cx+Au=0. 
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But this is still wrong, because it does not take account of the restriction that x must 
be nonnegative. So if any x; = 0, the trial solution can still be a minimum as long as the 
partial derivative is nonnegative. We therefore have the condition 


p+Cx+Au-ve=0, 


where x> 0, v > 0, u is unrestricted, and each component of v must vanish unless the corre- 
sponding component of x vanishes. 

Now, given a basic feasible solution of the original constraints, one can compute the 
corresponding values of p + C x and can write down the modified conditions 


p+Cx+Au-v+Ez=0, 


where E is a diagonal matrix whose diagonal elements are +1, chosen to make z > 0 if 
u = v = O and x has the value associated with the given basic feasible solution of the original 
constraints. 

Wolfe then uses the normal simplex method for linear programming to minimize 
a : 

>: Zi» subject to these n + m constraints, treating the u; as free variables and having an 
k=1 

additional rule to the effect that no v, (x;) may be removed from the nonbasic set unless the 
corresponding Xj (v;) is in the nonbasic set. He shows that this restriction on the choice of 
variable to be removed from the nonbasic set cannot prevent him from reaching an optimal 
solution with > . * 0, which solves the original problem. 

It does not seem possible at present to make a comprehensive summary of the relative 
advantages of Wolfe's and Beale's algorithms, but some things are clear. 

The chief advantage of Wolfe's algorithm is that very little modification of the simplex 
code for linear programming on a computer is needed to turn it into a code for quadratic 
programming. The SHARE linear-programming code for the IBM 704 required modification in 
some ten instructions for this purpose. Beale's algorithm would certainly need more than ten 
new instructions, though it should not prove a major programming problem, given a simplex 
code for linear programming. There is the further point that Beale's method does not seem 
suited to the inverse matrix method, for a variety of reasons, notably the fact that the number 
of basic variables plus the number of rows in the expression for C always exceeds the number 
of nonbasic variables. This remark does not of course apply to the variant proposed at the 
beginning of Section 3 when r<< n- m, 

The advantages of Beale's algorithm are: (1) that it is more direct and seems more 
natural; (2) it is more compact, particularly when the rank of the quadratic part of the objective 
function is small; (3) arising out of the first point is the fact that extensions of the algorithm are 
more easily developed. It is not clear how many of the problems considered in Section 3 can 
conveniently be handled by Wolfe's method. 

A vital point is the number of iterations required in the two methods. Wolfe's method 
may have an advantage here. Wolfe's method requires m iterations at the beginning to elimi- 
nate the u-variables from the nonbasic set, which correspond in Beale's method to the prelim- 
inary task of expressing C as a function of the nonbasic variables only. And at least n - m 
further iterations are needed to make all the Zi nonbasic. Some of these iterations may not be 
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needed in Beale's method if the first trial solution is a good one. But with Beale's method it is 
certain that if a new nonbasic restricted variable has to be introduced the final solution will not 
be obtained until C is again in standard form. In general this will involve removing all the 
present nonbasic free variables from the nonbasic set. Wolfe's method does not seem to suffer 
from this difficulty. I do not believe that this is a major handicap on the standard problem. 

But if one wishes to solve the parametric problem, with the linear part of the objective function 
multiplied by a parameter and the solution required for all values of the parameter, I believe 
that the "long form" of Wolfe's algorithm, described in [5], which is designed specifically for 
this problem, will prove superior unless it exceeds the maximum convenient storage capacity 
of one's computer. 
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An optimization problem which frequently arises in applications of 
mathematical programming is the following: 


F(A. B)= min [f) (x)) +f, (x2) +... +f), 
BL RES co SES 

where f; are convex functions. In this paper, the function F is studied 

and shown to satisfy F(A, B) = M(A) + N(B), where M and N are 

increasing and decreasing convex functions, respectively. Also, the 

functional equation F(A, C) = F(A, B) + F(B,C) - F(B,B) is estab- 

lished. These results generalize to the continuous case 


= 
F(A, B) = min { f(t, x (t)) dt, 


° 


with x(t) increasing and A <x(0) <x(T)<B. The results obtained in 
this paper are useful for reducing an optimization problem with many 
variables to one with fewer variables. 











1. INTRODUCTION 

A general type of problem not infrequently encountered in applications of mathematical 
programming is the following. It is desired to determine a projected sequence of activity 
levels X14, X9, +++, X,, Over a succession of time periods i= 1, 2, ..., m comprising a given 
planning interval 1<i<m. Associated with any allowable sequence of levels x, is a known 


1 
cost F(x;, Xo, e+e, x); furthermore, the function F is additively separable in the individual 





variables Xi, .é., 
m 
F(x;, Xo, soe, Xn) = f, (x) + fy (Xp) © ass +&) z z. f; (x;) . 
i=1 


Finally, the individual functions fi are convex. The problem then is to determine the activity 
m 

levels so as to minimize ze f; (x;) ‘ 
i=1 





*Manuscript received January 6, 1958. 
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This optimization is to be achieved relative to a set of linear restraints (equations or inequali- 
ties) which characterize the particular programming problem at hand; a sequence x;| is 
"allowable" (or "feasible") if it satisfies these restraints, implying that such a projected 
sequence of activity levels is physically realizable. 

If any sequence of values Xj is plotted as a function of the time i, then it appears as 
a series of increasing! and decreasing portions with corresponding peak and trough values. 
Now suppose that it is possible to argue from the mathematics of the problem,? or it is 
desirable to assume, that the peaks and troughs must occur at fixed times j, k, etc. That is, 


a feasible {x4} must satisfy, e.g., 





Xy < Xo eee 


lA 


Hi1 SX, 
* 2 Muy 2 eee 2 Mey 2 > etc. 


If this is so, the problem of minimizing F(x,, Xo, +++5 Xn) can be modified as follows. First, 
solve separately each of the problems 


G(A, B) = min fo(X5) +... + £,_ 4 (x, | 
ae) <p (2! 2) j-1 5-1 

2 j-1 
H(B, C) = min - isa a) teen + feo 2) : etc.; 


B =X 41> eee .%.9 = 
second, solve the problem 


(1) min ty (A) + G(A, B) + f9(B) + H(B, C) + | . 


with the variables A, B, C, ... (i.e., the renamed variables Xy> Xj, Xs + .) subject to 
alternating inequalities 


A<B, | a. ae etc. 


There is one apparent disadvantage to this reduction of the problem—while the original 
function F to be minimized is additively separable in X1, Xo, +++, the function that appears in 
(1) is not additively separable in its variables A, B,... . It is the main purpose of this paper 
to show that this difficulty is only apparent. It is, in fact, true that the function G can be 
written as a sum G, (A) + Go (B) with G,, Gp convex; similarly for H(B, C), etc. 

Let us state the result obtained in this paper more completely, with the notation to be 
employed. Let A, B be any real numbers with A < B. Let fi fo, ee ‘. be convex functions 
of a single real variable, each defined for all real values and bounded from below. Let 








linc reasing (decreasing) will always mean nondecreasing (nonincreasing). 
2For such a case see W. Karush and A. Vazsonyi, "Mathematical Programming and Employ- 
ment Scheduling,'' Nav. Res. Log. Quart. 4:297-320 (1957). 
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(2) F(A,B)= rea. 9 fn) + fg (xq) +... + £,(x,)]- 


Then it follows that? 


(3) F(A, B) = M(A) + N(B), ASB... 


where M is an increasing convex function and N is a decreasing convex function. A further 
useful result is the functional equation 


(4) F(A, C) = F(A, B) + F(B,C) - F(B, 8B), ARABS €.. 


In fact, we shall derive (3) from (4). 

There are analogous results for the continuous analog of the above discrete problem. 
In this case, the sequence of values Xj is replaced by a function x, with values x(t) defined 
for all t such that 0< t < T, and the sequence of functions f; is replaced by f(t, x) defined 
for OS t <T, -w< x < +0. Then 


T 


(5) F(A, B) = min f f(t, x@)) dt , with x(t) increasing 
O 


and A<x(0)<x(T)<B. 


With appropriate "smoothness" conditions on f(t, x), the results (3) and (4) remain valid in 
this continuous case. We shall content ourselves with a brief outline of the proof in this case 
in the final section of the paper. 


2. PRELIMINARIES 

It will be convenient to collect here some definitions and some known properties of 
convex functions. 

Let f be a real-valued function of a single variable x, defined for all real x. Then f 
will be taken to be convex in the following sense: for any x, y, Pp, q, 


(6) f(px + gy) < pf(x) + af(y), 
whenever p+q=1,p20,q20. 


Geometrically, this means that the line segment joining two points of the graph of f never falls 
below the graph. It is well known that the convexity of f implies the continuity of f; it also 
implies that the derivative 





Df(x) = lim it + 0) - $00 
t>0 t 





3The case of increasing x; also takes care of the case of decreasing levels, i.e., 


> i i i Ss « = =w = - = - 
C2 yy i eae zy, 2 D. Simply make the substitutions A D, x) Yuet Xy Yy> B G.. 
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exists, except possibly on a denumerable set of points. The function Df is increasing and 
consequently has left-hand and right-hand limits for every x which satisfy 








(7) lim_ Df(x +t) < lim is Df(x +t). 
t>0 t>0 


In addition, the right-hand derivative D*f and the left-hand derivative D°f exist everywhere; 


f(x + t) - f(x) D*t(x) = lim f(x + t) - f(x) 
t F t>0* t 








D' f(x) = lim _ 
t>0 


Both Df and D*f are increasing, and they are equal, respectively, to the left-hand and right- 
hand sides of (7). Where Df(x) exists, we have D f(x) = D*f(x) = Df(x). 

We say that g(x) is convex on P < x in case it is defined there and condition (6) holds 
for x, y in that range. In this case, we must restrict the properties stated above to the open 
half-line P < x; in particular, g may have a discontinuity at x = P. 

Let us consider now some specialized definitions convenient for our work. Let {x} be 
an increasing sequence, for example, Xj SX_ S.-- SX. Then the subrange r<i<s ofthe 
full range 1< i<n is said to be a maximal range* for Xj in case 


2X but x 


< < 
xX, and %, < His 


s’ r-1 


(the latter is to hold whenever 1 < r-1 or s+1 <n). We say that a sequence x; is (A, B) 
admissible in case A <x, < Xp <... <x, <B; thus, F(A, B) is the minimum of f; (x;) over 
(A, B) admissible sequences. This minimization problem we term Problem (A, B). Finally, 
a sequence fe has Property (A, B) in case, by definition, the sequence is (A, B) admissible 
and the following holds: let r <i<s be a maximal range on which the sequence has the 
constant value C; we require 


Ss 
if A<C<B, then >. D*£,(C) 2 0 
i=j 


j 
if A<C<B, then >. D°f,(C) < 0 
i=r 


for any j satisfying r <j <s. In essence, Property (A, B) guarantees that if we try to raise 
(lower) a right-hand (left-hand) section of {x;| lying over a maximal range, then we will not 


n 
decrease the cost function pam f; (x;) 5 
1 





For simplicity in the discussion to follow we write formulas as if r < s—the reader can make 


the modifications necessary for r = s. It may be, for example, that all maximal ranges are 
single values. 
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LEMMA 1: 


PROOF: For convenience, we shall use X, Y, Z to denote vectors of n components. 


3. THE FUNCTIONAL EQUATION 

In this section we shall establish the functional equation (4), following the proof of two 
preliminary lemmas. Recall that we are assuming that the functions fi are convex. In this 
section, however, we shall not require that they be bounded from below. 


A necessary and sufficient condition for 
that {%;| have Property (A, B). 
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Let X, Y be (A, B) admissible. Let 


XQ) = X+AZ, where Z=Y-X, 


Then X(A) is a family of (A, B) admissible vectors with 


Let 


where, of course, 


Now, a necessary and sufficient condition for X to be minimizing is that 


for all Z such that X + Z is (A, B) admissible. Since G(A) is convex in \—this is imme- 
diately derivable from the definition (6) and the convexity of the f; —the last requirement is 


equivalent to 


We have 


for all except possibly a denumerable set of A. Taking the limit as »A>07, and using the 
results of the preceding section, we obtain the following criterion for X to be minimizing: 


X(0) = X, X(1) = Y. 


GQ) = F(X()), 
n 

F(x) =)" f(x). 
i=1 


G(0) =~ = min G(A) 
eta 33 


0 < D*tG(0). 


n 
DG(A) = >. D f; (x; + dZ;) ° Zi 
i=] 










{i} to solve Problem (A, B) is 
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n 

(8) o< - D*£, (x) + 4, with D*(D") for z, > 0(<0), 
i=l 


for all Z such that X + Z is (A, B) admissible. 

Now, let us argue that if X solves Problem (A, B) then it has Property (A, B). Suppose 
X solves Problem (A, B), then condition (8) holds. Let r< i <s be a maximal range of X on 
which x has the constant value C. Let j satisfy r< j <s. Suppose, first, that C< B. Put 
























z= A, for j- ics, 
= 0, otherwise. 
For arbitrary positive A, sufficiently small, X + Z is (A, B) admissible. From (8) we obtain 
Ss 
ee ye D*i,(C) 20. 
i=j 


Next, suppose that A<C. Put 


z, = -A, for 2< i <j, 


This time (8) yields 


Thus, Property (A, B) is established. 
Now, suppose that X has Property (A, B). To establish condition (8) it will be con- 
venient to rewrite that condition. Let 


Zi> for Zz; > 0, 
Zz, = 

0, otherwise , 

Zi for 2; < 0, 
2; = 

Oo, otherwise . 


Then 
for all i. 
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+ 
Dj 


i 
D f; (x;) ’ 


S* (j, k) 


k 
ba D; , (all upper or all lower signs). 
i=j 


Then 


n n n 
+ " + + - .- 
¥ D* £,(x;) + 2 = > Di Zz; + > D; 2; - 
1 1 1 


Hence, to establish (8) it is sufficient to establish 
> s 
+ + - - 
D; 2; 2 0 and Ds D; Zz 2 0 
i=r i=r 


for any maximal range of X (since X can be decomposed into a series of constant sections 
over maximal ranges). A little manipulation presents these inequalities in the form 


s 
S*(r,s) + zt + pa S*(i,s) + (zy - zi 1) > 0 
i=r+1 
(9) 
s-l 
z S"(r,i) - (2; - 2541) + S (r,s) - Zs > 0. 
i=r 


The vector X + Z is (A, B) admissible by assumption, hence it is increasing. Since X 
has a constant value C over a maximal range, over this range we must have 


< €C cos e 
S, = Sng ¢ = 4, 
That is, 
+ - + - , 
s . < . : r<i<e-t. 
& * & = Sg * “es? =o 


A little consideration shows that since zy is nonnegative, Za, nonpositive, and at least one of 


Zr 2, is 0, the last inequality implies 


+ + - 
(10) te | and ZS 441 > 


i e<b<g@=t. 


Further, if C = A, then Zi = 0 in the maximal range; and if C = B, then Zz = 0 in the maximal 
range. 

We are now in a position to verify conditions (9). Suppose, first, that A < C< B. By the 
definition of Property (A, B) we see that each of the expressions S*(S") appearing in (9) is 
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nonnegative (nonpositive). By (10) it follows that inequalities (9) hold. Suppose, next, that 
C=A. The first inequality (9) follows from (10) and Property (A, B), while the second 
inequality follows from Zi = 0 onthe maximal range. Finally, the case C = B is argued 
analogously. This completes the proof of Lemma 1. 


The next lemma relates the minimizing problems for two different intervals (A, B). 
LEMMA 2. Suppose M<N and that Y solves Problem (M, N). Let 
MB<A<B <8. 


Then the solution X of Problem (A, B) is given by: 


m= A, ify, < A : 
= > if A<y,< B, 
x =B, if B< y, 


PROOF: By Lemma 1, Y has Property (M, N); by the same lemma it suffices to show 
that X, as defined above, has Property (A, B). Let r<i<s be a maximal range of X 
on which it takes the constant value C. Suppose, first, that A <C< B. Then Y shares 
this maximal range with X and has the value C on this range. Thus, Property (A, B) 
for X, in this case, follows from Property (M, N) for Y. Next suppose that A =C < B. 
(We may clearly assume A < B—the lemma is trivial for A= B.) Let 


Cy < Cg < ee < G 


be the constant values of Y not exceeding A with corresponding maximal ranges 


resis, | a Ae) oe 
Then % 3 i< Sy is the maximal range of X for the value A. It is required to show that 
Sp 
9 + P 
(12) >. D*4,(A) > 0, ry<j<8,- 
i=j 


Consider any fixed j, and choose the q for which 


Condition (12) becomes 


Sq41 


Sq Sp 
ta D*£,(A) + } a D*£,(A) +... + >: D*#,(A) > 
i=j i q+1 isr), 


Vv 
Oo 


1=r 
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By Property (M, N) for Y, and the monotonic increasing character of the functions pt fi, we 
obtain 


Sq Sq 
0 < b D*£(C,) < > D*#,(A) , 
i=j i=] 
Si Sk 
0< >: D*4(C,) < rb D*4,(A) , a . 
i=r), i=r), 


This establishes (12). The argument for the case A < C = B is similar. This completes the 
proof of Lemma 2. 

We are now in a position to establish the main result of this section. 

THEOREM 1: Let fj, i= 1, 2, ...,n, be n convex functions of a single real variable 
defined for all values of the variable.> Let F(A, B) be defined by (2). Then for arbitrary A, 
B, C satisfying A < B<C, F satisfies the functional equation (4). 


PROOF: Let X be minimizing for Problem (A, C). Let I, J, K be subsets of the index 
set {1, 2, ...,n} defined by 


I = all i with x; < B, 
J = all i with x; = B, 
K = all i with x, > B. 


Since X is minimizing, 


F(A, C) = th f; (x;) + bm f; (B) + a f; (x;) ‘ 
iin I iim J iin kK 
By Lemma 2, 
F(A, B) = Z f; (x;) * Pe f; (B) + ,> f; (B) 
i in I iin J iin K 
= i. f.(B (2) . 
eons i in I - * ee; i) + ee i) 





°The restriction that the f. be defined for all values of the real variable x is onlya convenience. 
If, for example, the functions are defined and convex on P <x, we can easily extend the 
definitions to all x preserving convexity. 
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By definition of F, As 
n 
F(B,B) = ) | £(B). 
i=1 
Equation (4) is now readily checked to complete the proof. Sun 
AN EXAMPLE 
We pause at this time for a simple example illustrating the functional equation (4). 
* 
Suppose each (convex) f; has a unique (finite) minimum point x;. Thus, f; is strictly 
7 * 
decreasing on -o <x< Xj and strictly increasing on x, < x <a. Suppose further that 
o * . 
Xj < Xp <-+e < K- 
For intervals (A, B) which contain all the x , the minimizing vector is clearly x* with 
components x; . Otherwise, let 
‘" < 4 < < " < < ‘ 
%j-1 <A x; wien X, B £ Xa: 
Then by Lemma 2 the minimizing vector is 
x, = A for i = 1,2,...,j-1, oes 
x =x for j<i<k | 
i i =? a ™ 9 
x =B for i=k+1, k+2,...,n. a 
This result is intuitively evident. (If x. < A, then x; = A; if B< x then xX; = B.) We have 
j-1 k n (14) 
os 
F(A, B) = rd £,(A) + >: f,(x;) + >: f, (B) . 
i=l i=j i=k+1 The 
To verify equation (4) we take B < C and define p, q by 
ar 
at <4, tu. «+ Cem. . 
Observe that 
" . The 
(13) either k=p-1 or k+1l=p-41, and X41 = %-1 = B- 
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As with F(A, B), we have 


p-1 q n 
F(B, C) = >: {,(B) + >: f,(x;) + ‘> f,(C) . 
i=1 i=p q+1 


Summing and using (13), we obtain 


j-1 k q n 
F(A, B) + F(B,C) = ‘ £,(A) + ‘io f,(x;) + >: f,(x;) + Zz f,(C) 
j p q+1 
n p-1 
+ >. {,(B) + pe f, (B) 
k+1 iz] 
j-1 q ‘ n n 
. >. i,(A) + 7. f,(x;) + >: f,(C) + Zz. f, (B) 
1 j q+1 i 


F(A, C) + F(B, B), 
as desired. 


5. DECOMPOSITION OF FUNCTION 


We return to the main argument, and deduce the decomposition of F, shown in (3), from 
the functional equation (4). 


LEMMA 3: Suppose that g(x) is continuous and convex for P < x. Let 


(14) G(P, x) = min g(y), Pes. 
es 9<z 


Then g(P, x) is continuous and convex in x for P < x, 


PROOF: The continuity of G(P, x) follows readily from the continuity of g(x). To 
argue convexity, let 


P< 2; < Zo, p+qe-=tl p>0, q>0. 
There exist X14, Xo such that 
g (x,) = G(P, Z1) > P £ xy = 21 
g (Xo) - G(P, Zo) ’ P< Xo = Zo + 





Recall that convexity on P <x implies continuity only on 
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From the convexity of g, 


pG(P, z,) + qG(P, Zz) = pg(x,) + ag(x9) > g(Px, + GX). 


But, P< PX; + GXq < PZ, + GZo; consequently, by the definition of G, 
g(px; + 4X5) > G(P, pz, + dz»). 
Hence, 
pG(P, z,) + qG(P, 25) > G(P, pz, + QZ), 


as required for the convexity of G(P, x). 
In preparation for the next lemma we define 


Fj (P, x) = - _— pee ose ¥ ie) » eR saw as 


In particular, 
F(A, B) = F(A; B). 


We see that 
F,(P,x) = _ min f, (y) 
P<y <x 
(15) 
F, P, x) = min F. (P, + f, ‘ 
jot (P» ») » fis | j (Py) ja 


LEMMA 4: For any P, F(P, x) is continuous and convex in x for P< x. 


PROOF: We make the proof by induction on the sequence F,. Fix P. By Lemma 3, 


F, (P, x) is continuous and convex for P < x. Now assume that F; ip, x) has the same property. 


It follows that 

F; (P, x) + fiat (x) 
also has this property. By (15) and Lemma 3 again, the property is preserved for Fyat (P, x). 
Since F = Fy this completes the proof. 


In the remainder of this section we assume the existence of a constant K such that? 


f;(x) > K, -0o < X < +0, oo eee 





7Notice that the assumption of boundedness from below does not imply that f; has a finite 
minimizing point X55 OB, f(x) = e*. 
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The function F may be extended to infinite values of the arguments as follows: 






















i} 


F(-«, x) lim F(A, x) 


A->-o 


F(x, +0) = lim F(x, B) 
B+ +0 


F(-w,+0) = lim F(A,B). 
A—>-o@ 
B+>+0 


For x fixed, F(A, x) is decreasing as A decreases and bounded from below. Hence, F(-«, x) 
is well defined for each x. In a similar manner it may be verified that F(x, +o) and F(-0, +00) 
are well defined. 


LEMMA 5: F(-o, x) (F(x, +00)) is continuous, monotonic decreasing (increasing) and 
convex in x for -« <x <+, Also, the equation (4) holds with A = -« or C = +0. 


PROOF: The validity of (4) for the stated infinite values follows by allowing the proper 
variables to tend to +o in (4). 


For each (finite) A, F(A, x) is monotonic decreasing and convex for A < x; the 
monotonicity follows from the definition of F, and the convexity from Lemma 4. Thus 


F(A, px + gy) < pF(A, x) + qF(A, y), A<x<y, 
with p+q=1, p>0, q2>0. Letting A—-wo, we obtain the convexity of F(-o, x). Continuity 


is a consequence of convexity, while monotonicity follows from the same property for F(A, x). 
To argue similar properties for F(x, +0) we may proceed as follows: 


‘Z “a 
F(x, B) = min | > f. (x.) | 
mt... 2a48(3 °% 3 
ag eee © J 
= min b g.(y,) | (where, g,(z) = £,_,(-2)) 
“B <y, <2. ER | 1 


G(-B, -x). 


From the preceding, G(-B, t) is a monotonic decreasing convex function of t for -B< t. Hence, 
G(-B, -x), i.e., F(x, B), is a monotonic increasing convex function of x for x <B. The prop- 
erties of F(x, +o now follow as above. 

We are now prepared to derive (3). By (4), and the second statement of Lemma 5, 


F (-00, 0) 


F(-o, B) + F(B, +0) - F(B, B) 


F(-o, B) + F(A, ow) - F(A, B). 
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Hence (3) holds with 








(The 
1 
M(A) = F(A, ©) - ; F (-«, «) 
N(B) = F(-«, B) - . F(-«, «) . 
THEOREM 2: Let the f; (x) be convex functions for -«© < x<« which are bounded 
from below. Let F(A, B) be defined by (2). Then F(A, B) may be written in the form 
(16) F(A, B) = M(A) + N(B), A<B 
where M(x) and N(x) are, respectively, increasing and decreasing convex functions for 
-o <x <o. Conversely, if F is any function of the form (16), then there are functions fi such 
that F is expressible in the form (2); explicitly, 
Ther 
F(A, B) = min | M(x,) + whee) ‘ 
A <x, £ Xp <B L 
PROOF: The first part of the statement follows from the paragraph preceding the 
theorem and Lemma 5. The converse holds because, by the monotonicity of M and N, 
min M (x;) + N (Xo) = M(A) + N(B) = F(A, B). 
A <% <z—< 8B 
1 2 
6. ALGORITHM FOR MINIMIZING POINT Hen 
The last section accomplishes our main goal for the discrete problem. Before passing 
on to the continuous analog we discuss a method of obtaining a minimizing vector X for the 
discrete problem, which is a special application of dynamic programming.® (The method whic 
actually does not require convexity of the f;-) Although the algorithm does not bear directly min: 
on the main goal, it will serve to answer a question that might arise in the reader's mind as to 
an effective procedure for the determination of the optimal vector. con: 
Suppose given continuous functions fy) fo, eas « on the domain A < x < B. Define — 
continuous functions Hj; F; on the domain as follows: con: 
H, (x) = f, (x) , F, (x) = min H, (y) line 
Ho(x) = F,(x) + f,(x) , F(x) = min Ho(y) 
2 1 2 2 A<y<x 2 THI 
in t 
H,-1 (*) = F.,-2 (%) + f,-1 (*) f F-1%) = ‘ min H,,_ 1 (%) 
HE @ = F434) + £6 F(x) = min BE tx). 
a) = Fog) +) 5 FG) =. min HAC 
8B ellman, R., ''The Theory of Dynamic Programming,'' Am. Math. Soc. Bull. 60:503-516 (1954). 
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(The present F; (x) is the earlier F; (A, x).) The required minimum value is given by 


F(A, B) = F,(B). 


Construct a sequence xy Xa-v Sissety Xy inductively as follows: 
x minimizes H,, (x) on A<x<B 

Xl minimizes Hy-1 (*) on A < x <i, 
xX minimizes H, (x) on A < X < Xp. 


There may be many such sequences; for any one we have, by the definition of H; and Fj 


F,(B) = Fy_1(%,) + £,(%,) 


Fy-1(%) = Fa-o(%-a) + fyi py) 


F,(X_) = £,(%,). 
Hence, 
F,, (B) = f, (X,) + fo (Xo) arene £(%) > 


which proves that (x), Xo, eeu x) is a minimizing point. This construction yields all 
minimizing points. 

The algorithm just described requires the construction of the functions H, and F.. This 
construction is graphically straightforward; the monotonic decreasing function F, is the "lower 
monotonic envelope" of H.. If the f, are convex, so are the H, and F., and the graphical 
construction of these functions, as well as of the values x., is simple. If the f, are piecewise 
linear, that is, have polygonal graphs, then the functions F,, H, are of the same type, and the 
algorithm is especially simple. (Here, convexity is not being assumed.) 


THE CONTINUOUS PROBLEM 
The proof of properties (3) and (4) in the case of the continuous problem (5) is outlined 
in this section. Consider 


min ft, x), t fixedin 0<t<T. 
-0O < xX < +0 
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Let x = m(t) solve this problem: 


f(t, m(t)) = min f(t, x), 
x 


We suppose that m(t) is a continuous function of t and consists of a finite number of increasing 
and decreasing sections on 0 <t <T. 

Now, we consider any (A, B) admissible x(t), that is, any increasing function on 
0<t <T with A <x(0), x(T) < B for which the integral in (5) is well defined. (We are 
assuming adequate regularity of the function f(t, x).) It can be argued that, without increasing 
the integral in (5), it is possible to replace any feasible x(t) by one with the following prop- 
erty: there is a partition of 0 < t < T into a finite number of subintervals such that on each of 
these "maximal" subintervals x(t) is either coincident with m(t) or constant. In considering 
the minimization problem (5), it is thereby sufficient to restrict attention to (A, B) admissible 
functions with this property. 

We may now define a Property (A, B) analogous to the property for the discrete case 
as follows: let r <t <s be a maximal subinterval on which x(t) = C; if C < B, then 


s 
il fr(t,C) dt > 0 
u x 


for every u satisfying r <u <s; if A<C, then 


u 
J f(t, A) dt < 0. 


In this definition fF (f,) designates the right-hand (left-hand) partial derivative of f with 
respect to x (which exists for every t, x by virtue of the fact that f is convex in x for each t). 
After this definition, the argument follows the strategy of proof for the discrete case. 
Lemma 1 holds as stated; Lemma 2 is also valid where the derived solution x(t) equals A 
whenever y(t) < A, equals B whenever y(t) > B, and equals y(t) otherwise. The functional 
equation (4) now follows from Lemma 2 as in Theorem 1. Next, the monotonicity and convexity 
of F(-o, x) and F(x, +) are established—the proof is somewhat different since we cannot 
use induction as in Lemma 4—and from this equation (3) is argued from equation (4) as before. 























PROBLEMS 

It has been suggested that this journal might serve as a medium for publishing problems 
in the area of Logistics — with the idea that other persons might be interested in those problems 
and might submit comments. Readers are invited to submit brief statements on applied and 
theoretical problems in Logistics. Address letters to Managing Editor, Naval Research Logis- 
tics Quarterly, Office of Naval Research, Washington 25, D. C. 


The following letter has been received: 


Dear Sir: 

This letter is in response to the recent letter by Ralph E. McShane and Henry Solomon. 

In their letter they offer certain criticisms of work which was partly done by Robert J. Aumann 
and myself. 

The criticisms by McShane and Solomon seem to me to be the following: 

(1) Wrong use of the MIP (Military Improvement Plan) list. 

(2) The preference relation might in some circumstances be intransitive (e.g., a pre- 
ferred to b preferred to c preferred to a). 

(3) One equipment will seldom be applicable to more than one system. 

(4) One or two minor remarks. 

Criticism (1) is clearly their major point, to which they devote about half their letter. 

When Aumann and I first became interested in this problem, one of the first questions we 
came up against was the intention and significance of the “MIP Priority list”. I wish to make the 
following points: 

(A) The meaning of the “priorities” given in the list is not made clear in any precise 
way in the list itself (which is promulgated by very high naval authority). 

(B) This list was not designed to be used in a numerical way with a mathematical 
scheme for solving assignment problems, but was designed for the guidance of human beings 
making decisions in the manner traditional in the Navy: namely, by using best judgment. 

(C) Though the word “priority” has a fairly specific meaning in mathematical economics, 
it is not clear that the high naval authorities who promulgated the list used the word in that 
sense. Rather, I believe, they used the word “priority” in a much more general sense, meaning 
that improvements with higher priority are more important, but not specifying how much more 
important. 








NOTES 


(D) Interpretation of the MIP list in the sense which McShane and Solomon seem to 
desire leads to assignment decisions which are obviously and greatly inferior to decisions 
which are reached by a more common sense interpretation. I refer here to real situations, 
not hypothetical ones. The reason for this is clear; the interpretation which McShane and 
Solomon suggest, the technical mathematical economics meaning, is that a higher priority is 
infinitely more important than the very next lower priority. 


(E) Aumann and I suggest finding out how much more important a higher priority is 





than a lower one by the same reasonable procedure which we suggested for determining other 
values; namely, by obtaining expert opinion on what to do in fairly simple though realistic 
situations, and finding out what values these decisions imply. Thus, if the MIP list really is 
intended to mean that each priority is infinitely more important than the next lower priority, 
then the expert opinion (which we actually obtained, though in a small way) would surely have 
reflected this, as the experts (presumably) have a real feeling for what the MIP list means, 
even though they do not express themselves in economic terms. Our results were not consis- 
tent with the infinite importance difference between priorities suggested by McShane and 
Solomon. 

With regard to criticism (2), of course we might run into intransitivity. Whenever you 
go up to a human being and ask “which do you prefer, a or b; which do you prefer, b or c; which 
do you prefer, c or a?” it is possible he will say “a better than b, b better than c, c better than 
a”. Surprisingly enough, in more complicated cases, this actually happens, though not too often. 


But, this risk exists as soon as you consult expert opinion on anything. The risk is no greater 





with our technique that any other which relies on expert opinion. Moreover, despite the remarks 
in the letter by McShane and Solomon, the risk is no greater when one equipment is applicable 
to many systems than when one equipment is applicable to one system. 

Criticism (3) is not considered to be valid. It often happens that one equipment is appli- 
cable to many systems, though perhaps with reduced effectiveness in some of them. This isa 
characteristic feature of the assignment problem which we were considering, and is one of the 
important factors which contributes to its complexity. 

One of the “minor remarks” listed under (4) is that the terminology used is complex and 
not too clear. I fully agree. It is complex because the problem is complex, and nothing can be 
done about that. We worked hard to simplify and clarify the terminology, but I admit that we 
did not perhaps, succeed too well in this direction. One of the difficulties is that certain terms 
are fixed by historical use, and may not be lightly changed. 

I am pleased to see that the work to which Aumann and I contributed has stirred up inter- 
est. It is of course clear that further work is needed. However, I believe that the work already 
done has considerable potential usefulness, needing only to be applied to fulfill its potential. 

Sincerely, 


JOSEPH B. KRUSKAL /sgn/ 
University of Michigan 

Ann Arbor, Michigan 

9 April 1959 
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RECENT PUBLICATIONS 


THEORIE DES GRAPHES ET SES APPLICATIONS, by Claude Berge. Dunod, Paris, 
France, 1958. viii + 275 pp. 

This is the second book ever written on graph theory. The only previous book is the 
classical D. Konig, “Theorie der endlichen und unendlichen Graphen,” 1936; reprinted by 
Chelsea Publications, Inc., 1950. Since that time there has been a considerable amount of 
significant research in this area, and a resurgence of interest in the theory of graphs. As the 
author notes in his introduction, graph theory has been recently applied to the behavioral 
sciences, information theory, cybernetics, game theory, transportation networks, and linear 
programming, as well as the theory of sets, matrices, and other purely abstract disciplines. 
Thus, this book is a welcome arrival on the mathematical scene, for it contains a wealth of 
useful, recent, and interesting material. 

The publishers have done a beautiful printing job and the author uses figures skillfully 
in his exposition. Unfortunately, however, the book appears to be directed at professional mathe- 
maticians, and is not likely to be readily understood by behavioral scientists without considerable 
effort. The logical organization of the material is unorthodox; for example, the basic graph 
theoretic concepts of distance and tree are not introduced until Chapters 12 and 16 respectively 
(out of 21 chapters). 

Without any doubt, the references constitute the weakest part of the book. Among other 
things, there are theorems attributed to the wrong person, and theorems named after incongru- 
ously associated pairs of authors. The entire Appendix V, which was written by J. Riguet, is 
more or less obtained by translating the powerful enumeration theorem of Pélya and one of his 
examples from German into French; nevertheless, Pdlya is only credited with one definition. 

The first three appendices are on linear programming and game theory, and are original 
with the author. Appendix I indicates the interaction between game theory and graph theory. 
Appendix II, on transportation problems, makes strong use of purely graph theoretic results in 
its application to some well known problems in linear programming, including the assignment 
problem and the transportation problem. Appendix III, on the utilization of the notion of poten- 
tial for transportation networks, is an extension of the preceding appendix and discusses the 
traveling salesman problem. Finally, Appendix IV consists of a good collection of 14 unsolved 
problems in graph theory. Altogether, the five appendices constitute an important portion of the 
book since they provide definite evidence of the importance of graph theory to other disciplines. 
The author gives frequent examples from several applied fields, and his original researches 
which are included in the book are of a high caliber. 


(Reviewed by F. Harary) 


RECENT PUBLICATIONS 


THE ARMY AND ECONOMIC MOBILIZATION, by R. Elberton Smith. United States 
Government Printing Office, 1959. 728 pp., index, $5.25. 

This volume, one of the series UNITED STATES ARMY IN WORLD WAR I, is the 5th 
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During World War II the War Department became one of the principal agencies of the 
Government in administering, as well as planning, the nation’s economic mobilization. This 
volume presents a comprehensive view of the impact of war on the national economy and gives 


the reader a fund of knowledge that will enable him to understand the complex tasks associated 


with Army procurement and economic mobilization in World War II. It should be of particular 


interest to those in the fields of economics and/or logistics. 
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